Про техники оптимизации

Поучительная история о техниках оптимизации наглядно.

Техзадание

Объявим в рамках топика небольшой конкурс по архитектурно-ориентированной оптимизации программного обеспечения.
Вкратце, код состоит из пачки функций, производящих невнятные на первый взгляд манипуляции с исходными данными, и примочки-драйвера, который n раз запускает неоптимизированную версию, затем оптимизированную, сравнивает насчитанные циферки, и, в случае их совпадения, выдает отношение времени выполнения. Вот так:

Executing original code… done
Executing optimized code… done
Checking results… PASSED
Number of runs: 3
Original code average time: 11.954537 sec.
Optimized code average time: 1.052994 sec.
Speedup: 11.35

Разрешено использовать любые техники оптимизации, компилятор GCC с любыми опциями, и, скажем, сервер с двумя четырехъядерными процессорами Intel Xeon E5420 2.5 GHz.
Вот, кстати, код:

/*
* contestopt.c:
*/

#include "contest.h"

void function3(double a[][NUM], double b[][NUM], double c[][NUM]);
void function2(double a[][NUM]);
double function4(unsigned int second);

double ***p;
double ****x;
double ****y;
double ****z;
double ***wrk1;
double ***wrk2;
double ***bnd;
static double a[NUM][NUM], b[NUM][NUM];

void function1(double result[][NUM], double gos[2], unsigned int second)
{
  function2(a);
  function2(b);

  function3(a, b, result);
  gos[0] = function4(second);
}

void function3(double a[][NUM], double b[][NUM], double result[][NUM])
{
  int i, j, k;

  for (i = 0; i < NUM; i++) {
    for (j = 0; j < NUM; j++) {
      for (k = 0; k < NUM; k++) {
        result[i][j] = result[i][j] + a[i][k] * b[k][j];
      }
    }
  }
}

void function2(double a[][NUM])
{
  int i, j;
  double first = 0.0;
  double second = 1.0;
  a[0][0] = first;
  a[0][1] = second;
  for (i = 0; i < NUM; i++) {
    for (j = 0; j < NUM; j++) {
      if (!(i == 0 && (j == 0 || j == 1))) {
        a[i][j] = first + second;
        first = second;
        second = a[i][j];
      }
      if (j % 20 == 0 && j != 0) {
        first = first * (j + 1) / (NUM);
        second = second * (j + 1) / (NUM);
      }
    }
    first = ((i + 1) * NUM) / first;
    second = ((i + 1) * NUM) / second;
  }
}

void function5(double ***A, double val)
{
  int i, j, k;
  for (i = 0; i < XX; i++)
    for (j = 0; j < YY; j++)
      for (k = 0; k < ZZ; k++)
        A[i][j][k] = val;
}

void function6(unsigned int second)
{
  int mimax = XX;
  int mjmax = YY;
  int mkmax = ZZ;

  int imax = mimax - 1;
  int jmax = mjmax - 1;
  int kmax = mkmax - 1;

  int i, j, k, l;
  double val;
  srand((unsigned int) second);
  p = (double ***) malloc(sizeof(double **) * XX);
  wrk1 = (double ***) malloc(sizeof(double **) * XX);
  wrk2 = (double ***) malloc(sizeof(double **) * XX);
  bnd = (double ***) malloc(sizeof(double **) * XX);

  x = (double ****) malloc(sizeof(double ***) * XX);
  y = (double ****) malloc(sizeof(double ***) * XX);
  z = (double ****) malloc(sizeof(double ***) * XX);

  for (i = 0; i < XX; i++) {
    p[i] = (double **) malloc(sizeof(double *) * YY);
    wrk1[i] = (double **) malloc(sizeof(double *) * YY);
    wrk2[i] = (double **) malloc(sizeof(double *) * YY);
    bnd[i] = (double **) malloc(sizeof(double *) * YY);

    x[i] = (double ***) malloc(sizeof(double **) * YY);
    y[i] = (double ***) malloc(sizeof(double **) * YY);
    z[i] = (double ***) malloc(sizeof(double **) * YY);

    for (j = 0; j < YY; j++) {
      p[i][j] = (double *) malloc(sizeof(double) * ZZ);
      wrk1[i][j] = (double *) malloc(sizeof(double) * ZZ);
      wrk2[i][j] = (double *) malloc(sizeof(double) * ZZ);
      bnd[i][j] = (double *) malloc(sizeof(double) * ZZ);

      x[i][j] = (double **) malloc(sizeof(double *) * ZZ);
      y[i][j] = (double **) malloc(sizeof(double *) * ZZ);
      z[i][j] = (double **) malloc(sizeof(double *) * ZZ);

      for (k = 0; k < ZZ; k++) {
        x[i][j][k] = (double *) malloc(sizeof(double) * 4);
        y[i][j][k] = (double *) malloc(sizeof(double) * 3);
        z[i][j][k] = (double *) malloc(sizeof(double) * 3);

      }

    }
  }

  val = (double) rand() / (10.0 * RAND_MAX);
  for (i = 0; i < XX; i++)
    for (j = 0; j < YY; j++)
      for (k = 0; k < ZZ; k++) {
        for (l = 0; l < 3; l++) {
          x[i][j][k][l] = 1.0;
          y[i][j][k][l] = val;
          z[i][j][k][l] = 1.0;
        }

        x[i][j][k][3] = 1.0;
      }
  function5(bnd, 1.0);
  function5(wrk1, 0.0);
  function5(wrk2, 0.0);
  for (i = 0; i < imax; i++)
    for (j = 0; j < jmax; j++)
      for (k = 0; k < kmax; k++)
        p[i][j][k] = (double) (k * k) / (double) (kmax * kmax);
}

void function7()
{
  int i, j, k;
  for (i = 0; i < XX; i++) {
    for (j = 0; j < YY; j++) {
      free(p[i][j]);
      free(wrk1[i][j]);
      free(wrk2[i][j]);
      free(bnd[i][j]);

      for (k = 0; k < ZZ; k++) {
        free(x[i][j][k]);
        free(y[i][j][k]);
        free(z[i][j][k]);
      }
      free(x[i][j]);
      free(y[i][j]);
      free(z[i][j]);
    }
    free(p[i]);
    free(wrk1[i]);
    free(wrk2[i]);
    free(bnd[i]);

    free(x[i]);
    free(y[i]);
    free(z[i]);

  }
  free(x);
  free(y);
  free(z);

  free(p);
  free(wrk1);
  free(wrk2);
  free(bnd);
}

double function4(unsigned int second)
{
  int nn = NN;
  double gosa, s0, ss;
  int mimax = XX;
  int mjmax = YY;
  int mkmax = ZZ;
  int imax = mimax - 1;
  int jmax = mjmax - 1;
  int kmax = mkmax - 1;
  int iCount, jCount, kCount, loop, i, j, k;

  function6(second);
  for (loop = 0; loop < nn; loop++) {
    gosa = 0.0;
    for (k = 1; k < kmax - 1; k++)
      for (j = 1; j < jmax - 1; j++)
        for (i = 1; i < imax - 1; i++) {
          s0 = x[i][j][k][0] * p[i + 1][j][k]
            + x[i][j][k][1] * p[i][j][k]
            + x[i][j][k][2] * p[i][j][k]
            + y[i][j][k][0] * (p[i + 1][j + 1][k] - p[i + 1][j - 1][k]
                      - p[i - 1][j + 1][k] + p[i - 1][j - 1][k])
            + y[i][j][k][1] * (p[i][j + 1][k + 1] - p[i][j - 1][k + 1]
                      - p[i][j + 1][k - 1] + p[i][j - 1][k - 1])
            + y[i][j][k][2] * (p[i + 1][j][k + 1] - p[i - 1][j][k + 1]
                      - p[i + 1][j][k - 1] + p[i - 1][j][k - 1])
            + z[i][j][k][0] * p[i - 1][j][k]
            + z[i][j][k][1] * p[i][j - 1][k]
            + z[i][j][k][2] * p[i][j][k - 1] + wrk1[i][j][k];
          ss = (s0 * x[i][j][k][3] - p[i][j][k]) * bnd[i][j][k];
          gosa = gosa + ss;
          wrk2[i][j][k] = p[i][j][k] + 0.1 * ss;
        }
    for (iCount = 1; iCount < imax - 1; iCount++)
      for (jCount = 1; jCount < jmax - 1; jCount++)
        for (kCount = 1; kCount < kmax - 1; kCount++)
          p[iCount][jCount][kCount] = wrk2[iCount][jCount][kCount];
  }
  function7();
  return gosa;
}

double exponent(int i, double x)
{
  double z = 1.0;
  int j;
  for (j = 0; j < i; j++) {
    z = z * x;
  }
  return z;
}

double function8(float *a, float p)
{
  int i, j;
  double h, x;
  double s = p;
  h = 1.0 / VAL;
  for (i = 1; i <= VAL; i++) {
    x = h * (i - 0.5);
    for (j = 0; j < ORDER; j++) {
      s = (a[j] * exponent(j, x)) + s;
    }
  }
  return s * h;
}


* This source code was highlighted with Source Code Highlighter.

Тестируется на значениях:
#define NUM 600
#define XX 65
#define YY 33
#define ZZ 33
#define NN 50
#define VAL 10000000
#define ORDER 30
#define TOLERANCE 1e-8

Примерный вид драйвера:
/*
* driver.c: Contest driver.
*/

#include <sys/time.h>
#include <stdio.h>
#include <stdlib.h>
#include <inttypes.h>
#include <time.h>

#include "contest.h"
#include "hpctimer.h"

double _RES[NUM][NUM] = { 0.0 };
double _GOS[2] = { 0.0, 0.0 };
double RES[NUM][NUM] = { 0.0 };
double GOS[2] = { 0.0, 0.0 };
float A[ORDER] = { 0.0 };

double function8(float *a, float p);
void function1(double result[][NUM], double gos[2], unsigned int second);

int is_results_eq(double _res[][NUM], double _gos[2], float _p,
         double res[][NUM], double gos[2], float p)
{
  int i, j;

  if (fabsf(_p - p) > 1E-1) {
    printf("_p p\n");
    printf("%f\n", _p);
    printf("%f\n", p);
    return 0;
  }

  if (fabs(_gos[0] - gos[0]) > 1E-1) {
    printf("_gos gos\n");
    printf("%f\n", _gos[0]);
    printf("%f\n", gos[0]);
    return 0;
  }

  for (i = 0; i < NUM; i++)
    for (j = 0; j < NUM; j++) {
      if (fabs(_RES[i][j] - RES[i][j]) > TOLERANCE) {
        printf("[%d, %d] =\n%f\n%f\n", i, j, _RES[i][j], RES[i][j]);
        return 0;
      }
    }
  return 1;
}

int main(int argc, char **argv)
{
  unsigned int second;
  unsigned int i, j;
  float p = 9.0;
  double _res, res;
  int nruns = 3;
  hpctimer_time_t t0, t1;
  double torig = 0.0, topt = 0.0;
  
  srand(time(NULL));
  second = 1219304613 + (1000 - (2000.0 * (double) rand() / (double) (RAND_MAX + 1.0)));
  for (i = 0; i < ORDER; i++) {
    A[i] = (float)i;
  }

  for (i = 0; i < NUM; i++) {
    for (j = 0; j < NUM; j++) {
      _RES[i][j] = 0.0;
      RES[i][j] = 0.0;
    }
  }

  /* Original code */
  printf("Executing original code ... ");
  fflush(stdout);
  t0 = hpctimer_gettime();
  for (i = 0; i < nruns; i++) {
    _function1(_RES, _GOS, second);
    _res = _function8(A, p);
  }
  t1 = hpctimer_gettime();
  torig = hpctimer_getdiff(t0, t1) / nruns;
  printf("done\n");
  fflush(stdout);

  /* Optimized code */
  printf("Executing optimized code ... ");
  fflush(stdout);
  t0 = hpctimer_gettime();
  for (i = 0; i < nruns; i++) {
    function1(RES, GOS, second);
    res = function8(A, p);
  }
  t1 = hpctimer_gettime();
  topt = hpctimer_getdiff(t0, t1) / nruns;
  printf("done\n");
  fflush(stdout);

  printf("Checking results ... ");
  fflush(stdout);

  if (!is_results_eq(_RES, _GOS, _res, RES, GOS, res)) {
    printf("RESULTS ARE DIFFERENT!\n");
  } else {
    printf("PASSED\n");
    printf("Number of runs: %d\n", nruns);
    printf("Original code average time: %.6f sec.\n", torig);
    printf("Optimized code average time: %.6f sec.\n", topt);
    printf("Speedup: %.2f\n", torig / topt);
  }
  return 0;
}


* This source code was highlighted with Source Code Highlighter.

Ускорение, Правило Первое

Об оптимизации на уровне исходного кода можно было бы утверждать: «Не важно, что делает код — важно, как он это делает». Но выведенное мною Первое Правило Оптимизации несколько противоречит этому утверждению:

Первое Правило
Курить матчасть и алгоритмы

Ибо при первом взгляде на код становится ясно, что некоторые простые действия он выполняет через анальное отверстие.
Так вот, даже беглого взгляда хватает, чтобы с омерзением выпилить функцию exponent, которая значительно садит количество полезных телодвижений (запустив какой-нибудь профилировщик, можно обнаружить, что она жрет 80% общего времени счета). Возведение в степень, которым она занимается, в function8 легко реализуется динамически.
В самой function8 подсознательно просматривается вычисление интеграла, но просторы для ее оптимизации обсудим позже.
Разностная схема, используемая в function4, не сильно-то оптимизируется переписыванием кода, но, вынеся накопление gosa в отдельный цикл, все же можно ее ускорить.
Умножение матриц в function3 разгоняется множеством различных способов, но на первом этапе достаточно заметить, что перемножаемые матрицы абсолютно идентичны, и сократить расходуемую в этой функции память в два раза. А то и вовсе нагуглить изящный алгоритм возведения матрицы в квадрат. Я уже не говорю про Штрассена, Копперсмита-Винограда и прочих извращенцев.

Ускорение, Правило Второе

Я не знаю, специально ли автор кода напихал в него столько избыточности, но основная часть полученного ускорения обязана Второму Правилу:

Второе Правило
Уничтожить структуры данных, не несущие смысловой нагрузки

Правило, кстати, прекрасное, и работает не только на оптимизацию, но и на удобочитаемость, и на лучше-компилируемость. Лирически цитируя одну замечательную женщину, «нет смысла использовать рекурсию там, где ее нет смысла использовать».
В этом смысле function4 нагружена прямо-таки до неприличия. С ужасом заметив, чем инициализируются четырехмерные массивы x, y и z, уничтожить их немедленно. За ними могут отправиться wrk1 и bnd. Я уже не говорю о таких мелочах, как mimax, iCount, ss и прочих излишествах. В итоге код становится куда короче и эстетичнее.
function6 и function7 целиком ликвидируются, остатки полезных действий не нагрузят function4 ни логически, ни по времени, поэтому безболезненно переносятся в нее. function5 внезапно самоуничтожается. В function3, как уже говорилось, можно использовать одну исходную матрицу, а function2 вызвать единожды.
Уже гораздо проще дышать, а получаемый speedup начинает вдохновлять на подвиги.

Ускорение, Правило Третье

Третье Правило
Минимизировать ошибки предсказания переходов

Сейчас все объясню. Это техника оптимизации ветвлений. Скомпилированный цикл или оператор ветвления выглядит примерно так (сильно упрощенно):
1 mov eax, ebx
2 jne $a
3 mov ecx, ebx
4 jmp b
5 a: ...
6 b: ...

Представим, что происходит в это время в конвейере. Инструкции обрабатываются в несколько этапов (Fetch — Decode — Execute — Write-back), и после загрузки инструкции 2 должен работать модуль предсказания переходов, потому что мы пока еще не знаем результат этой операции и не можем выбрать следующую.
Что происходит, когда модуль предсказания переходов «ошибается»? Правильно — конвейер разгружается. А это накладные расходы. Уменьшить количество ошибок модуля предсказаний переходов мы можем, уменьшив количество обращений к нему.
Как?
  • размещать наиболее вероятные ветви в начале ветвлений
  • выносить выше (по уровню вложенности в циклах) инвариантные ветвления
  • использовать разворачивание циклов
Поясню последнее: если вместо
for (i = 0; i < n; i++) {
  sum += p[i];
}

написать
for (i = 0; i < n; i+=4) {
  sum += p[i];
  sum += p[i + 1];
  sum += p[i + 2];
  sum += p[i + 3];
}

мы уменьшим количество предсказаний переходов в четыре раза.

Есть у этой техники пачка подводных камней, например, кратность количества итераций разворачиваемым кускам, оптимальная длина разворачивания и прочие. Но в каждом индивидуальном случае можно найти баланс и отхватить лишние проценты ускорения этим относительно простым способом.
Пример разворачивания цикла можно увидеть в оптимизированной версии function3.
Еще неплохо бы добавить вынесение из циклов инвариантных ветвлений, которые очень мозолят глаза в function2, но не так-то просто их оттуда выцарапать. Кто реализует лучше, чем я — тот молодец.

На почитать: Р. Гербер, А. Билл «Оптимизация ПО: сборник рецептов». 2010, Intel

Ускорение, Правило Четвертое

Четвертое Правило
Эффективно работать с памятью

В частности, с кэш-памятью. Есть такое свойство локальности ссылок — свойство программ повторно/часто обращаться к одним и тем же объектам; и имеет место временная (temporal) локализация и пространственная (spatial) локализация. Кэширование идет с упреждением, и на этом основании есть смысл обращаться по последовательным адресам, чтобы избежать промахов по кэшу.
В том же перемножении матриц, строка располагается в памяти последовательно, а столбец рассеивается, и известный эффект ускорения дает транспонирование второй матрицы.
Раз уж речь зашла, приведу еще пример: если структура (struct) полностью не помещается в кэш, и при этом содержит внутри себя массив, лучше использовать массив структур — вместо массива в структуре.
А еще, чтобы написать переносимую программу, не стоит привязываться к размеру кэша.
И многие другие тонкости. Почитать — Крис Касперски «Техника оптимизации программ. Эффективное использование памяти». 2003, БХВ-Петербург

Ускорение, Правило Пятое

Какая у нас там архитектура процессора?

Пятое Правило
Векторизовать вычисления, где это возможно

Самое время использовать SSE, благо, векторизовать можно практически любой код, работющий с массивами.
Небольшой ликбез.
Архитектура SSE (x86):
  1. MMX — обработка целочисленных векторов (64-битные регистры)
  2. SSE — работа с вещественными и целочисленными векторами (128 бит)
Инструкции могут работать как со скалярами, так и с упакованными векторами. Например:
     mulss xmm1, xmm0  |  mulps xmm1, xmm0
XMM0 4.0 3.0 2.0 1.0   |  4.0  3.0  2.0  1.0
                  *    |   *    *    *    *
XMM1 5.0 5.0 5.0 5.0   |  5.0  5.0  5.0  5.0
                  =    |   =    =    =    =
XMM1 4.0 3.0 2.0 5.0   |  20.0 15.0 10.0 5.0

Инструкции SSE умеют копирование, арифметику, сравнение, поразрядные операции и еще всякие хитрые преобразования.

Как их использовать? Есть ассемблерные векторные инструкции, есть специальные классы C++, есть, в конце концов, автоматическая векторизация — самая переносимая. Мне захотелось использовать интринсики.
Intrinsics — набор функций и типов данных, поддерживаемых компилятором, для предоставления высокоуровневого доступа к SSE инструкциям.
Подключаем xmmintrin.h, пишем векторизованный код, компилируем с флагом -msse3. Компилятор распределяет регистры XMM#, планирует команды и режимы адресации. Удобно, интуитивно просто.
Например, сложить одной инструкцией четыре числа типа float:
#include <xmmintrin.h>
void add(float *a, float *b, float *c)
{
  __mm128 t0, t1;
  t0 = _mm_load_ps(a);
  t1 = _mm_load_ps(b);
  t0 = _mm_add_ps(t0, t1);
  _mm_store_ps(c, t0);
}

В данном конкретном случае переписывается на интринсики перемножение матриц (заодно оно становится блочным). function8 в принципе тоже векторизуется, но тут вообще особый случай, так что пока поиспытываю ваше терпение.
На этом этапе код начинает сниться, а каждый вырезанный из толстого цикла оператор дает ощутимый скачок производительности.

Ускорение, Правило Шестое

Смотреть на граничные условия. Заменять структуры данных на более простые, где это возможно. Сводить вычисления к задачам, для решения которых существуют быстрые алгоритмы.

Шестое Правило
Использовать моск

В драйвере результаты перемножения матриц сравниваются с точностью до TOLERANCE, а остальные два числа — с точностью до десятых. Методом научного тыка выясняется, что в function8 для достижения требуемой точности достаточно две трети интервала, а первые три миллиона (всего-то) итераций можно отбросить. Мало того, достаточно учитывать каждую шести-семитысячную итерацию.
Получаем
double function8(float *a, float p)
{
  int i, j, val, step;
  double h, x, z;
  double s = p;

  h = 1.0 / VAL;
  val = 0.33 * VAL;
  step = val / 500;
  for (i = val; i <= VAL; i += step) {
    x = h * i;
    z = 1.0;
    for (j = 1; j < ORDER; j++) {
      z *= x;
      s += j * z;
    }
  }
  return s * h * step;
}

и радуемся. Нравится? Мне — нет. Ибо можно вообще избавиться от операции возведения в степень.
double function8(float *a, float p)
{
  unsigned int j;
  double s = p / VAL;
  for (j = 0; j < ORDER; j++) {
    s += (a[j] / (j + 1));
  }
  return s;
}

Получается сложность O(ORDER) — по сравнению с исходной O(VAL) * O(ORDER^2). Того же результата можно достичь, использовав числа Бернулли (которые уже не так тривиально, но все-таки можно заметить, выписав алгоритм вручную).
Можно обратить внимание на особенности заполнения возводимой в квадрат матрицы и увидеть нули в четных строках.
Можно заменить трехмерные массивы на одномерные, адресуя их наподобие:
m1 = (YY — 1) * (ZZ — 1); m2 = (ZZ — 1); p[i * m1 + j * m2 + k]. Эта небольшая хитрость позволяет использовать memset для заполнения многомерных матриц.
Можно и даже нужно использовать __inline для небольших функций.

И еще. Поскольку в данном случае имеется в распоряжении восьмиядерная машина, нельзя не использовать OpenMP. На распараллеливание function3, function4 и function8 даже много фантазии не потребуется. Чуть менее очевидна возможность использования параллельных секций в function1 (ибо function3 и function4 не зависят по данным).

На сладкое — версия оптимизированной программы. На идеал не претендует, но ускорение в честных несколько сотен раз выдает.
/*
* contestopt.c:
*/

#include <sys/time.h>
#include <omp.h>
#include "contest.h"

void function3(double a[][NUM], double c[][NUM]);
void function2(double a[][NUM]);
double function4(unsigned int second);

static double p[XX][YY][ZZ] __attribute__ ((aligned(64)));
static double wrk2[XX][YY][ZZ] __attribute__ ((aligned(64)));
static double a[NUM][NUM] __attribute__ ((aligned(64)));

void function1(double result[][NUM], double gos[2], unsigned int second)
{
  #pragma omp parallel sections
  {
    #pragma omp section
    function3(a,result);

    #pragma omp section
    gos[0] = function4(second);
  }
}

void function3(double a[][NUM], double result[][NUM])
{
  int i, i2, j, j2, k, k2;
  double *restrict rres;
  double *restrict rmul1;
  double *restrict rmul2;
  unsigned int BS = 8;
  int remainder = NUM % 8;
  int limit = NUM - remainder;

  function2(a);
  if (!(NUM & 1)) {
    #pragma omp for private(k, j, i, i2, j2, k2, rres, rmul1, rmul2)
    for (i = 0; i < limit; i += BS)
      for (j = 0; j < limit; j += BS)
        for (k = 0; k < limit; k += BS)
          for (i2 = 0, rres = &result[i][j], rmul1 = &a[i][k]; i2 < BS;
                         ++i2, rres += NUM, rmul1 += NUM) {
            _mm_prefetch (&rmul1[8], _MM_HINT_NTA);
            if (!(NUM & 1))
              for (k2 = 0, rmul2 = &a[k][j]; k2 < BS; ++k2, rmul2 += NUM) {
                __m128d m1d = _mm_load_sd (&rmul1[k2]);
                m1d = _mm_unpacklo_pd (m1d, m1d);
                j2 = 0;
                __m128d m2 = _mm_load_pd (&rmul2[j2]);
                __m128d r2 = _mm_load_pd (&rres[j2]);
                _mm_store_pd (&rres[j2],
                _mm_add_pd (_mm_mul_pd (m2, m1d), r2));
                j2 +=2;
                m2 = _mm_load_pd (&rmul2[j2]);
                r2 = _mm_load_pd (&rres[j2]);
                _mm_store_pd (&rres[j2],
                _mm_add_pd (_mm_mul_pd (m2, m1d), r2));
                j2 +=2;
                m2 = _mm_load_pd (&rmul2[j2]);
                r2 = _mm_load_pd (&rres[j2]);
                _mm_store_pd (&rres[j2],
                _mm_add_pd (_mm_mul_pd (m2, m1d), r2));
                j2 +=2;
                m2 = _mm_load_pd (&rmul2[j2]);
                r2 = _mm_load_pd (&rres[j2]);
                _mm_store_pd (&rres[j2],
                _mm_add_pd (_mm_mul_pd (m2, m1d), r2));
              }
          }
  } else {
    #pragma omp for private(k, j, i, i2, j2, k2, rres, rmul1, rmul2)
    for (i = 0; i < limit; i += BS)
      for (j = 0; j < limit; j += BS)
        for (k = 0; k < limit; k += BS)
          for (i2 = 0, rres = &result[i][j], rmul1 = &a[i][k]; i2 < BS;
                        ++i2, rres += NUM, rmul1 += NUM) {
            _mm_prefetch (&rmul1[8], _MM_HINT_NTA);
            for (k2 = 0, rmul2 = &a[k][j]; k2 < BS; ++k2, rmul2 += NUM) {
              for(j2 = 0; j2 < BS; j2++ ) {
                rres[j2] += rmul1[k2] * rmul2[j2];
              }
            }
          }
  }
  if (remainder) {
    #pragma omp for private(i,j,k)
    for (i = 0; i < limit; ++i)
      for (k = NUM - remainder; k < NUM; ++k)
        for (j = 0; j < limit; ++j)
          result[i][j] += a[i][k] * a[k][j];
    #pragma omp for private(i,j,k)
    for (i = limit; i < NUM; ++i)
      for (k = 0; k < NUM; ++k)
        for (j = 0; j < NUM; ++j)
          result[i][j] += a[i][k] * a[k][j];
    #pragma omp for private(i,j,k)
      for (i = 0; i < limit; ++i)
        for (k = 0; k < NUM; ++k)
          for (j = limit; j < NUM; ++j)
            result[i][j] += a[i][k] * a[k][j];
  }
}

void function2(double a[][NUM])
{
  int i, j ;
  double first = 0.0;
  double second = 1.0;

  __assume_aligned(a, 64);

  a[0][0]=first;
  a[0][1]=second;

  for (j = 2; j < NUM; j++) {
    a[0][j] = first + second;
    first = second;
    second = a[0][j];
    if( j%20 == 0 && j != 0) {
      first = first * (j + 1) / (NUM);
      second = second * (j + 1) / (NUM);
    }
  }
  first = NUM / first;
  second = NUM / second;
  for (i = 1; i < NUM; i++) {
    for (j = 0; j < NUM; j++) {
      a[i][j] = first + second;
      first = second;
      second = a[i][j];
      if( j%20 == 0 && j != 0 ) {
        first = first * (j + 1) / (NUM);
        second = second * (j + 1) / (NUM);
      }
    }
    first = ((i + 1) * NUM) / first;
    second = ((i + 1) * NUM) / second;
  }
}

__inline void function6(unsigned int second)
{
  int imax = XX - 1;
  int jmax = YY - 1;
  int kmax = ZZ - 1;
  int i, j, k;

  double kmaxDoubled = kmax * kmax;
  for (k = 0; k < kmax; k++)
    p[0][0][k] = (k * k) / kmaxDoubled;
  for (i = 0; i < imax; i++)
    for (j = 0; j < jmax; j++)
      for (k = 1; k < kmax; k++)
        p[i][j][k] = p[0][0][k];
}

double function4(unsigned int second)
{
  int nn = NN;
  double gosa = 0.0, val;
  int imax = XX - 2;
  int jmax = YY - 2;
  int kmax = ZZ - 2;
  int loop, i, j, k;

  function6(second);
  srand((unsigned int)second);
  val = (double) rand() / (10.0 * RAND_MAX);
  for (loop = 1; loop < nn; loop++) {
    for (k = 1; k < kmax; k++)
      for (j = 1; j < jmax; j++)
        for (i = 1; i < imax; i++) {
          wrk2[i][j][k] = 0.1 * ( p[i + 1][j][k] + p[i][j][k]
            + val * (p[i + 1][j + 1][k] - p[i + 1][j - 1][k]
                 - p[i - 1][j + 1][k] + p[i - 1][j - 1][k]
                 + p[i][j + 1][k + 1] - p[i][j - 1][k + 1]
                 - p[i][j + 1][k - 1] + p[i][j - 1][k - 1]
                 + p[i + 1][j][k + 1] - p[i - 1][j][k + 1]
                 - p[i + 1][j][k - 1] + p[i - 1][j][k - 1])
            + p[i - 1][j][k] + p[i][j - 1][k] + p[i][j][k - 1]);
        }
     for (i = 1; i < imax; i++)
      for (j = 1; j < jmax; j++)
        for (k = 1; k < kmax; k++) p[i][j][k] += wrk2[i][j][k];
  }
  for (k = 1; k < kmax; k++)
    for (j = 1; j < jmax; j++)
      for (i = 1; i < imax; i++) {
        gosa += p[i + 1][j][k] + p[i][j][k]
          + val * (p[i + 1][j + 1][k] - p[i + 1][j - 1][k]
               - p[i - 1][j + 1][k] + p[i - 1][j - 1][k]
               + p[i][j + 1][k + 1] - p[i][j - 1][k + 1]
               - p[i][j + 1][k - 1] + p[i][j - 1][k - 1]
               + p[i + 1][j][k + 1] - p[i - 1][j][k + 1]
               - p[i + 1][j][k - 1] + p[i - 1][j][k - 1])
          + p[i - 1][j][k] + p[i][j - 1][k] + p[i][j][k - 1];
      }
  return gosa;
}

double function8(float *a, float p)
{
  unsigned int j;
  double s = 0.0;

#pragma omp parallel
{
  #pragma omp for private(j) reduction (+:s)
  for (j = 0; j < ORDER; j++) {
     s += (a[j] /(j + 1));
    }
}
  return (s + p / VAL);
}


* This source code was highlighted with Source Code Highlighter.

Вот такой веселый и весьма полезный практический пример. Обобщим схему оптимизации кода:
  1. Алгоритмическая оптимизация
  2. Архитектурно-независимая оптимизация:
    • ввод-вывод в параллельный поток
    • распараллеливание
    • эффективная работа с памятью, эффективное кэширование
  3. Архитектурно-ориентированная оптимизация:
    • векторные инструкции
    • минимизация количества ошибок предсказания переходов
    • библиотечные функции от производителей процессоров
    • вынесение расчетов на видеокарту
    • etc.
И чтобы было все еще интересно. Приведенная программа driver.c содержит уязвимость, используя которую, можно получить ускорение в десятки тысяч раз, особо не напрягаясь. Хотя это уже будет эксплоит.

Спасибо за внимание.
Share post
AdBlock has stolen the banner, but banners are not teeth — they will be back

More
Ads

Comments 66

    +1
    Ого
      +13
      чо я себя идиотом почувствовал с IQ ниже, чем у креветки…
        +14
        Первая мысль — за такие имена функций и переменных неплохо было бы расстреливать…
          +10
          Да-да. По коду видно, что это какие-то вычисления. Автору было бы неплохо написать, что этот код решает такую-то задачу. Конечно можно по функциям сидеть и разбирать и найти известные алгоритмы, но у меня не возникает желания напрягать мозг для этого.
            +1
            первая, и единственная.
              +5
              в том-то и беда, что код мне достался именно в таком виде. и это нередкая ситуация из тех, когда приходится оптимизировать чей-то код
                +3
                Хорошо, что вы сохранили авторский стиль в названии функций. =)
                  +8
                  Я так полагаю, что автор уже мертв? :)
                  0
                  такие имена функций еще более менее читабельны. после обфускаторов/zend в еще большем лесе приходится копаться.
                    0
                    Вам придётся расстрелять ~всех выпускников мехмата МГУ. Подозреваю, и прочих матмехов. Иосиф Виссарионович, перелогиньтесь!
                      0
                      А можно разъяснить при чем тут мехмат МГУ? Их что-ли учат программировать в таком стиле?
                        0
                        Скорее, не учат программировать в другом. Подобный стиль написания кода похож на запись математических уравнений, которым учат матмехов, а вот умению писать читаемый код их не учит никто.

                        Я не то чтобы готов, на самом деле, искать причины; я просто неоднократно видел программный код в исполнении матмеха, и каждый раз он был похож на вышеприведённый, и даже хуже. Вот, красивый примерчик нашёл (разбивка на строки моя, форматирование оригинала):

                        miss+=fabs(X[i+(N+1)*(j+1)]+X[i+(N+1)*(j-1)]+X[i+1+
                        N+1)*j]+X[i-1+(N+1)*j]-4*X[i+(N+1)*j]-b[i+(N+1)*j]); 
                          0
                          Мда… А у нас вот наоборот, некоторые преподы лабы не принимали, если не было единого стиля форматирования кода и адекватных имен идентификаторов…
                            0
                            Зато наверняка и математика была слабее. Каждому — своё.

                            Ряд профессоров до сих пор на Фортране пишет. У программиста на Фортране своё, аутентичное понимание адекватности имён идентификаторов. В частности, длина имени должна быть не более 6 символов. Сложно при этом выдвигать промышленные требования к читаемости кода, и слава Богу, что они не пытаются этим заниматься.
                              +1
                              Насколько я знаю, более-менее современные версии фортрана уже таких требований не выдвигают… Или я не прав?
                            0
                            Как матмех (физтех) скажу, что сложно только до выделения X[i+(N+1)*(j+1)] в семантическую единицу. Потом формула становится очень хорошо читаемой. Как будто фокус на микроскопе настроил, и все стало хорошо видно.
                            Я так же (хотя это и выглядит как будто наоборот) с SQL поступаю — использую всегда полные имена «базаданных».«таблица».«поле» Текста получается много, зато всегда ясно, что откуда берется. А так как писать и переписывать SQL мне приходится много, то экономится много времени на том, что не надо смотреть, откуда же берется это поле. При чтении такого SQL тоже достаточно «настроить фокус» и все хорошо видно.
                      0
                      Отлично написано и наглядно. В статью в букмарки, книги в 2read.

                      p.s. нынче в большинстве случаев дешевле купить мощнее сервер, чем нанять толкового программиста который знает толк в оптимизации. или разрешить ему потратить только на это 1-5 дней.
                      p.p.s. еще есть злая штука — преждевременная и избыточная оптимизация. когда в момент разработки кода, когда еще толком не ясно что и как код будет делать — его пробуют улучшить и ускорить. в результате 2 дня оптимизируют функцию, которая потом выпиливается, т.к. agile и планы поменялись. (пример из жизни)
                        0
                        А можно узнать, что можно 2 дня оптимизировать в функции? Если Вы знаете оптимизирующие хаки языка, в котором пишите, то Вы либо сразу пишите код с оптимизацией, либо в течение 5 минут заменяете все что надо на оптимизированный вариант.
                        Пример с python
                        писать вместо str += 'vasya' str = ''.join([str,'vasya'])
                        Или вместо for i in range(5) i = 0; while i < 5: i+=1
                          +2
                          Можно, например функция прорисовки обьекта или какой-то кастомный контрол. Во время разработки черновика приложения допустимо использовать такое решение, которое займет меньше времени. Иначе, если оптимизировать все что пишешь — уйдет больше времени на придумывание хороших решений и к дню дедлайна будет сделана треть приложения. Отлично сделана — но за нее не заплатят.
                        +2
                        Все по делу, любой человек, который занимался оптимизацией подтвердит :)
                        Единственное замечание — попытка повторного использования счетчиков в циклах вкупе for private из open mp. Зачем, если можно объявлять счетчик в каждом цикле и он автоматически будет объявлен приватным при распараллеливании? Ну и заодно дать счетчикам осмысленные имена.
                          +1
                          скажем так, почти все оптимизации банальны, кроме SSE: все таки не каждый программист знает такие низкоуровневые функции, которые в большинстве приложений собственно и не требуются.
                          А данный код на том же фортране смотрелся бы куда изящнее.
                            0
                            Знаете, я бы не стал говорить обо всех оптимизациях.
                            Я бы даже сказал, что алгоритмические оптимизации зачастую очень сложно понять без специальных пояснений. В качестве примера, почитайте Майкла Абраша «Дзен графического программирования»: там очень наглядно показано, как из одного алгоритма получить совершенно другой потём серии оптимизаций.
                              0
                              Имелись в виду оптимизации рассмотренные в статье. А банальны, они в том смысле, что для тех кто этим занимается постоянно, приходится делать одни и те же оптимизации ( имею в виду научные вычисления).
                                0
                                в этом смысле статья больше предназначалась для тех, кто вообще оптимизацией не занимался. а так конечно, стоило последний абзац оставить и все.
                            +5
                            Меня сейчас пристрелят, но не могу не процитировать (произносится с дебильным выражением лица):
                            — «Папа, а с кем это ты сейчас разговаривал?» :)

                            По сути обсуждаемого вопроса: синтетические тесты тем плохи, что напоминают по полезности перекладывание кучки щебня в другую кучку детским совочком. В реальной жизни такая штуковины либо будет исполняться «как есть», пусть даже на бейсике и будет считать 30 секунд, а не 0.3, да и фиг с ней; если же надо считать часто и много — то начиная с подсчетов на GPU и заканчивая решениями на ПЛИС (в том числе и с использованием lookup tables), ибо прирост скорости будет несравнимый.
                              +3
                              Насколько я понимаю, то с чем занимался любовными утехами автор — это реальный проект. В котором надо было улучшить производительность. Видимо кого-то-таки не устраивало, что считалось 30 секунд вместо 0,3. И этот кто-то — начальство.
                                –1
                                Да тут вроде довольно-таки явно видно, что пример — тестовый. Я не верю в то, что будут писать ТАК реальную систему, там же непонятно нифига начиная от имен функций.

                                Ну, или просто у меня мало опыта, но я реальной жизни такого не видел еще. :)
                                  +1
                                  Всякое бывает. Помню года 4 назад пришлось мне разбираться в хтмплн-ом отчете с джаваскриптами, где все функции принимали параметр c именем O и назывались f1, f2, f3… Что при отсутствии типизации весьма доставляло.
                                  0
                                  А, да, еще: из общения с математиками у меня возникло ощущение, что пока они свои формулы не упростят донельзя, не оптимизируют, не обсудят с коллегами, и не посидят вместе еще пару дней, споря над результатом — на реализацию не выдадут. :)
                                    0
                                    А как иначе? Критерий успешности математика — красота теории, т.е. ее простота.
                              • UFO just landed and posted this here
                                  –2
                                  Всё это здорово, только при чем здесь высокая производительность?
                                  Из всего здесь написанного к делу имеет отношение только выбор правильных алгоритмов, который может повысить производительность в разы (точнее, на шаги сложности функции O(*) ). Неправильный выбор алгоритма — это баг, один из самых серьезных типов багов вообще. Этот баг не позволяет использовать программу на больших объемах входных данных (то есть в реальных условиях), и править этот баг нужно _обязательно_.
                                  Да, повышение локальности _данных_ (а не ссылок) тоже может в разы ускорить программу, но это опять же — правильный выбор структур данных и работающих с ними алгоритмов.
                                  Всё остальное… нет, ну не то чтобы неверно. Всё это правильно. Но просто эффект… Ну, допустим, вы, просидев пару месяцев с профайлером и книжкой Гербера-Билла, добились ускорения работы программы в целом… на 20%. А нагрузка на систему за эти пару месяцев выросла на 50%. Вы потратили кучу усилий, но копали не в ту сторону.
                                  «Низкоуровневая» оптимизация нужна и полезна, только не для высокой производительности, а для систем реального времени. Вот библиотеки Intel для работы, скажем, с изображениями — они просто бесценны, если ваша программа большую часть времени обрабатывает изображения. Какие-нибудь трансляции, или распознавание образов в реальном времени, или системы слежения за чем-то — и всё это на относительно слабом железе — вот в этих случаях вы будете выжимать из железа любые проценты. Но это «высокая удельная производительность».
                                  А в высоконагруженных системах по мере необходимости просто добавляется один-два-десять-сто новых серверов за 1000-1500 долларов каждый. И задача программиста — обеспечить безпроблемное горизонтальное масштабирование системы на эти новые сервера. Была обработка одной структуры данных на 20 серверах, стала на 25 — и всё это прозрачно, нигде ничего не останавливалось и не переиндексировалось. И вот там уже начинаются свои трюки, засады и просторы для улучшений.
                                    +2
                                    > добились ускорения работы программы в целом… на 20%

                                    Примерно в 200 раз.
                                    Посоветуете более подходящий блог?
                                      –4
                                      Почему не over9000? =)
                                      Я говорю о реальных проектах без явных ляпов, написанных нормальными программистами.
                                      Понятно, что можно придумать примеры, поддающиеся ускорению и в 200, и в 200 тысяч раз.
                                      Насчет блога — прямого попадания вроде нет, но по идее ближе всего блоги о Си++
                                      Реальность такова, что единственный смысл использовать сегодня Си или Си++ — как раз таки выжать максимум удельной производительности.
                                    +3
                                    самый отрезвляющий код за этот год.
                                      +6
                                      весьма лестно, учитывая, что пока еще только третье января.)
                                      0
                                      Последнего абзаца хватило бы, без простыни предшествующего текста.
                                        0
                                        1) А ты теперь скомпиль оба исходника в ICC с полной оптимизацией и сравни результаты.
                                        2) А как ты замеряешь время? По времени выполнения кода или по процессорному времени которое даётся потоку. К тому как я вижу ты юзаешь и системные функции, по этому тут уже очень много факторов будут. тебе надо будет считать процессорное время для kernel mode и для user mode. Потому что замеры обычного времени тебе ничего не дадут. Тем более для такого малого времени тестирования.

                                          0
                                            +1
                                            Подсчет по тактам процессора — это жесть. Потому что тут ты всё равно не учитываешь многозадачность. за 1 секунду, смена потока произойдет огромное кол-во раз и как на крути у тебя будут данные завышены и сильно привязаны к состоянию системы. Потому как процессоры новые имею в себе функции энергосбережения и автоматического уменьшения таковой частоты в зависимости от нагрузки. И по этому первичная калибровка тебе мало что даст. И опять же malloc и free — как не крути но тут будут вызваны функции системы по выделению памяти, а скорость этой операция напрямую зависит от загруженности памяти, т.к. будут затрачено время на поиск свободной страницы виртуальной и физической страницы (конечно это время сверх мало, но всёже вносит свои ограничения) и таким образом по немного но собирается куча обстоятельств которые вносят погрешности
                                              +1
                                              > Подсчет по тактам процессора — это жесть.

                                              Может и жесть, но точность высокая.

                                              С многозадачностью можно бороться запуская на чистой машине, где нет других «тяжелых» задач и путем многократного замера.

                                              Энергосбережение и уменьшение тактовой частоты не при делах, если в proc cpuinfo указан флаг constant_tsc (гарантирующий стабильность счетчика тактов). Ибо счетчик тактов (не только он) часто используется ОС-ми для реализации функций выдачи текущего времени. Разве что с turbo boost будут небольшие проблемы.

                                              Не каждый malloc/free вызывает ОС, более того, выделение-освобождение часто можно вынести за пределы измеряемого кода.
                                          0
                                          > Лирически цитируя одну замечательную женщину, «нет смысла использовать рекурсию там, где ее нет смысла использовать».

                                          Цитируя другого замечательного человека, «нет смысла использовать мутабельные переменные там где нет смысла их использовать».

                                          Кто-то не прав.
                                            +2
                                            знакомый код какой-то… такое ощущение, что его давали в качества задания на одной из «Зимних Школ Intel по параллельным вычислениям»…
                                              0
                                              А теперь вот это всё — взять и подсветить.
                                                0
                                                fixed
                                                0
                                                спасибо, для меня очень полезная статья
                                                  0
                                                  Код таймера дали выше. А можете залить куда-нибудь весь код в итоговом виде архивом? В драйвере используются функции _function1 и пр., откуда они?
                                                    +1
                                                      0
                                                      Мда, на ускорение посмотреть не получилось… Проц Core 2 Duo E8500, ось 32-бит. С включённым openmp неоптимизированный вариант выполнился нормально, оптимизированный застрял сам в себе. После отключения openmp оба варианта выполнились, но дали разные результаты. Тестировалось в gcc 4.4.3 и ~4.5.0.
                                                        0
                                                        занятно. перемножение матриц упало) попробуйте его упростить
                                                          +1
                                                          Предполагаю, что неоптимизированная версия считала на fpu (387), который немного точнее sse2. Лечить опциями "-msse2 -mfpmath=sse". На моей машине помогло.
                                                      0
                                                      #include <stdio.h>

                                                      int main() {
                                                      int n,i;
                                                      n = 33;
                                                      for (i = 0; i < n; i+=4) {
                                                      printf("%i %i %i %i\n", i,i+1,i+2,i+3);
                                                      }
                                                      return 0;
                                                      }

                                                      Выдаст не ожидаемый результат (он будет считать до 35), хотя я может не прав и в случае с array результат будет нульданные в C
                                                        0
                                                        чего-то я не понял, и почему результат неожидаемый? все верно… до i=32 (последний шаг), и выведет оно 32, 33 (32+1), 34 (32+2), 35 (32+3)
                                                          0
                                                          ну что получить надо до 32 (нет по коду ясно, что выдаст до 35) а получаем до 35. вопрос больше автору топика
                                                          0
                                                          хм, ну никто и не будет разворачивать цикл с нечетным количеством итераций по 4. упомянутый мною подводный камень с кратностью. берут потом остаток от деления и досчитывают…
                                                            0
                                                            хм, да возможно так лучше будет работать.
                                                            будет ли подобное работать при чтение файла, например?
                                                              0
                                                              я чет недопоняла. что должно работать? уменьшается количество обращений к модулю branch prediction — по факту уменьшения итераций в цикле.
                                                                0
                                                                Если не использовать системных вызовов (вход в который дороже ошибки предсказания в разы), и работать либо с хорошо буферизованным файлом либо с mmaped файлом, то будет работать.

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

                                                                branch prediction от раскрутки выигрывает мало, т.к. предсказание переходов просто обязано правильно предсказывать циклы. Неверных предсказаний на каждый заход в цикл будет от 0 (ибо переходы назад — предсказываются как происходящие) до 2-3 (после нескольких итераций предсказатель увидит, что ошибался и поверит в этот цикл). Ну еще +1 ошибка на выход из цикла. )

                                                                Что еще выигрывает от раскрутки — логика цикла (проверка пора ли выходить) срабатывает реже (не на каждую итерацию, а на итерации/ степень раскрутки)

                                                                Ссылки:
                                                                en.wikipedia.org/wiki/Loop_unwinding (общее + ссылки)

                                                                software.intel.com/sites/products/documentation/studio/composer/en-us/2009/compiler_c/optaps/common/optaps_hlo_unrl.htm (что говорит документация от компилятора Интел о полезности unroll)

                                                                software.intel.com/sites/products/documentation/studio/composer/en-us/2009/compiler_c/optaps/common/optaps_perf_optstrats.htm (общее из той же документации про циклы и их оптимизацию)
                                                            +1
                                                            насколько я помню, это вообще был итоговый конкурс, после зимней школы оптимизации. дали такой странный код, со словами: у вас есть 20 минут на его оптимизацию, юзайте всякие профайлеры и все что хотите. суть в том была, чтобы не сами алгоритмы оптимизировани на высоком уровне, а пользуюясь сторонними утилитами находить узкие места программы и переписывать их нормальным образом. насколько я помню, 80% кода вообще ничего не делала (что-то типа гоняла данные туда-сюда), достаточно было этот факт подменить, выкинуть эту часть и все — победа в кармане :) но за 20 минут реально сложновато в этом сраче разобраться.
                                                              –2
                                                              Почём сейчас коробок спичек?
                                                                0
                                                                они из неоптимизированной версии. да пожалуйста, тут
                                                                  +2
                                                                  да чтоб. промахнулась
                                                                  +2
                                                                  По поводу еще более серьезной оптимизации ветвлений и прочего есть вводная статья от IBM: http://www.ibm.com/developerworks/linux/library/l-gcc-hacks/

                                                                  Там правда с упором на GCC, ядро и промышленность, но тем больше доверия заслуживают методы
                                                                    +1
                                                                    Ну и серьезные дяди (и тети?) еще проводят в конце итеративный адаптивный тюнинг a-la acovea
                                                                      0
                                                                      пасиба, почитаем
                                                                      0
                                                                      все очень интересно,
                                                                      книгу Криса Касперски держал в руках, но так и не купил
                                                                      куплюобязательно

                                                                      еще раз спасибо за содержательную историю

                                                                      Only users with full accounts can post comments. Log in, please.