Оптимизация циклов C

Question

Я новичок в C из многих лет Matlab для численного программирования. Я разработал программу для решения большой системы дифференциальных уравнений, но я почти уверен, что сделал что-то глупое, так как после профилирования кода я был удивлен, увидев три цикла, которые занимают ~ 90% вычислений время, несмотря на то, что они выполняют самые тривиальные шаги программы.

Мой вопрос состоит из трех частей, основанных на этих дорогих циклах:

• Инициализация массива в ноль. Когда J объявляется как двойной массив, значения массива инициализируются нулем? Если нет, есть ли быстрый способ установить все элементы на ноль?

  void spam(){
      double J[151][151];    
      /* Other relevant variables declared */
      calcJac(data,J,y);
      /* Use J */
  }
  
  static void calcJac(UserData data, double J[151][151],N_Vector y)
  {
      /* The first expensive loop */
      int iter, jter;
      for (iter=0; iter<151; iter++) {
          for (jter = 0; jter<151; jter++) {
              J[iter][jter] = 0;
          }
      }
     /* More code to populate J from data and y that runs very quickly */
  }
  

• В процессе решения мне нужно решить матричные уравнения, определяемые P = I - гамма * J. Построение P занимает больше времени, чем решение системы уравнений, которую оно определяет, поэтому то, что я делаю, вероятно, ошибочно. В приведенном ниже относительно медленном цикле доступ к матрице, содержащейся в структуре «данные», является медленным компонентом или это что-то еще из цикла?

  for (iter = 1; iter<151; iter++) {
      for(jter = 1; jter<151; jter++){
          P[iter-1][jter-1] = - gamma*(data->J[iter][jter]);
      }
  }
  

• Есть ли лучшая практика для умножения матриц? В приведенном ниже цикле Ith (v, iter) - это макрос для получения iter-го компонента вектора, содержащегося в структуре N_Vector 'v' (тип данных, используемый решателями Sundials). В частности, есть ли лучший способ получить скалярное произведение между v и строками J?

  Jv_scratch = 0;
  int iter, jter;
  for (iter=1; iter<151; iter++) {
      for (jter=1; jter<151; jter++) {
          Jv_scratch += J[iter][jter]*Ith(v,jter);
      }
      Ith(Jv,iter) = Jv_scratch;
      Jv_scratch = 0;
  }
  

Answer 1

1) Нет, они не могут запоминать массив следующим образом:

memset( J, 0, sizeof( double ) * 151 * 151 );

или вы можете использовать инициализатор массива:

double J[151][151] = { 0.0 };

2) Ну, вы используете довольно сложный расчет для расчета положения P и положения J.

Вы вполне можете получить лучшую производительность. переходя по указателям:

for (iter = 1; iter<151; iter++) 
{
    double* pP = (P - 1) + (151 * iter);
    double* pJ = data->J + (151 * iter);

    for(jter = 1; jter<151; jter++, pP++, pJ++ )
    {
         *pP = - gamma * *pJ;
    }
}

Таким образом вы перемещаете различные вычисления индекса массива за пределы цикла.

3) Лучшей практикой является попытка вывести из цикла как можно больше вычислений. Как и в предыдущем цикле.

Answer 2

Во-первых, я бы посоветовал вам разделить ваш вопрос на три отдельных вопроса. Трудно ответить на все три; Я, например, мало работал с численным анализом, поэтому отвечу только на первый.

Во-первых, переменные в стеке не инициализированы для вас. Но есть более быстрые способы их инициализации. В вашем случае я бы посоветовал использовать memset:

static void calcJac(UserData data, double J[151][151],N_Vector y)
{
   memset((void*)J, 0, sizeof(double) * 151 * 151);
   /* More code to populate J from data and y that runs very quickly */
}

memset - это быстрая библиотечная подпрограмма для заполнения области памяти определенным шаблоном байтов. Просто так получилось, что установка всех байтов double в ноль устанавливает double в ноль, поэтому воспользуйтесь преимуществами быстрых подпрограмм вашей библиотеки (которые, вероятно, будут написаны на ассемблере, чтобы использовать преимущества таких вещей, как SSE).

Answer 3

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

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

Answer 4

Инициализация массива в ноль. Когда J объявляется двойным массив - это значения массива инициализируется до нуля? Если нет, есть ли быстрый способ установить все элементы в ноль?

Это зависит от того, где расположен массив. Если он объявлен в области видимости файла или как статический, то стандарт C гарантирует, что все элементы установлены в ноль. То же самое гарантируется, если вы установите первый элемент в значение при инициализации, то есть:

double J[151][151] = {0}; /* set first element to zero */

Устанавливая первый элемент в что-либо, стандарт C гарантирует, что все остальные элементы в массиве будут установлены в ноль, как если бы массив был статически размещен.

Практически для этого конкретного случая я очень сомневаюсь, что будет разумно выделить 151 * 151 * sizeof (double) байтов в стеке независимо от того, какую систему вы используете. Скорее всего, вам придется распределять его динамически, и тогда ничего из вышеперечисленного не имеет значения. Затем вы должны использовать memset (), чтобы установить все байты в ноль.

В Относительно медленный цикл ниже доступ к матрице, которая содержится в структуре «данных» медленный компонент или это что-то еще о петле?

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

Конечно, вы могли бы запутать код с помощью ручных оптимизационных операций, таких как обратный отсчет до нуля, а не вверх, или использование ++ i, а не i ++ и т. Д. И т. Д. Но компилятор действительно должен уметь обрабатывать такие вещи за вас.

Что касается сложения матриц, я не знаю математически наиболее эффективного способа, но подозреваю, что он не имеет никакого отношения к эффективности кода. Большой вор здесь - двойной тип. Если вам не нужна высокая точность, я бы рассмотрел использование float или int для ускорения алгоритма.