Решение систем линейных уравнений

8.2.2. Определение подзадач

При внимательном рассмотрении метода Гаусса можно заметить, что все вычисления сводятся к однотипным вычислительным операциям над строками матрицы коэффициентов системы линейных уравнений. Как результат, в основу параллельной реализации алгоритма Гаусса может быть положен принцип распараллеливания по данным. В качестве базовой подзадачи можно принять тогда все вычисления, связанные с обработкой одной строки матрицы A и соответствующего элемента вектора b.

8.2.3. Выделение информационных зависимостей

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

Для выполнения прямого хода метода Гаусса необходимо осуществить (n-1) итерацию по исключению неизвестных для преобразования матрицы коэффициентов A к верхнему треугольному виду.

Выполнение итерации i, 0<=i<n-1, прямого хода метода Гаусса включает ряд последовательных действий. Прежде всего, в самом начале итерации необходимо выбрать ведущую строку, которая при использовании метода главных элементов определяется поиском строки с наибольшим по абсолютной величине значением среди элементов столбца i, соответствующего исключаемой переменной xi. Поскольку строки матрицы A распределены по подзадачам, для поиска максимального значения подзадачи с номерами k, k<i, должны обменяться своими элементами при исключаемой переменной xi. После сбора всех необходимых данных в каждой подзадаче может быть определено, какая из подзадач содержит ведущую строку и какое значение является ведущим элементом.

Далее для продолжения вычислений ведущая подзадача должна разослать свою строку матрицы A и соответствующий элемент вектора b всем остальным подзадачам с номерами k, k<i. Получив ведущую строку, подзадачи выполняют вычитание строк, обеспечивая тем самым исключение соответствующей неизвестной xi.

При выполнении обратного хода метода Гаусса подзадачи выполняют необходимые вычисления для нахождения значения неизвестных. Как только какая-либо подзадача i, 0<=i<n-1, определяет значение своей переменной xi, это значение должно быть разослано всем подзадачам с номерами k, k<i. Далее подзадачи подставляют полученное значение новой неизвестной и выполняют корректировку значений для элементов вектора b.

8.2.4. Масштабирование и распределение подзадач по процессорам

Выделенные базовые подзадачи характеризуются одинаковой вычислительной трудоемкостью и сбалансированным объемом передаваемых данных. В случае когда размер матрицы, описывающей систему линейных уравнений, оказывается большим, чем число доступных процессоров (т.е. p<n ), базовые подзадачи можно укрупнить, объединив в рамках одной подзадачи несколько строк матрицы. Однако применение использованной в лекциях 6 и 7 последовательной схемы разделения данных для параллельного решения систем линейных уравнений приведет к неравномерной вычислительной нагрузке процессоров: по мере исключения (на прямом ходе) или определения (на обратном ходе) неизвестных в методе Гаусса для все большей части процессоров все необходимые вычисления будут завершены и процессоры окажутся простаивающими. Возможное решение проблемы балансировки вычислений может состоять в использовании ленточной циклической схемы (см. лекцию 6) для распределения данных между укрупненными подзадачами. В этом случае матрица A делится на наборы (полосы) строк вида (см. рис. 8.2):$$A=(A_0, A_2, \ldots, A_{p-1})^T, A_i = (a_{i_0}, a_{i_1}, \ldots, a_{i_{k-1}}), \; i_j=i+jp, \; 0 \le j < k, \; k=n/p$$

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

Распределение подзадач между процессорами должно учитывать характер выполняемых в методе Гаусса коммуникационных операций. Основным видом информационного взаимодействия подзадач является операция передачи данных от одного процессора всем процессорам вычислительной системы. Как результат, для эффективной реализации требуемых информационных взаимодействий между базовыми подзадачами топология сети передачи данных должны иметь структуру гиперкуба или полного графа.

8.2.5. Анализ эффективности

Оценим трудоемкость рассмотренного параллельного варианта метода Гаусса. Пусть, как и ранее, n есть порядок решаемой системы линейных уравнений, а p, p<n, обозначает число используемых процессоров. Тем самым, матрица коэффициентов А имеет размер nxn, и, соответственно, n/p есть размер полосы матрицы А на каждом процессоре.

Прежде всего, несложно показать, что общее время выполнения последовательного варианта метода Гаусса составляет:

Определим теперь сложность параллельного варианта метода Гаусса. При выполнении прямого хода алгоритма на каждой итерации для выбора ведущей строки каждый процессор должен осуществить выбор максимального значения в столбце с исключаемой неизвестной в пределах своей полосы. Начальный размер полос на процессорах равен n/p, но по мере исключения неизвестных количество строк в полосах для обработки постепенно сокращается. Текущий размер полос приближенно можно оценить как (n-i)/p, где i, 0<=i<n-1, есть номер выполняемой итерации прямого хода метода Гаусса. Далее после сбора полученных максимальных значений, определения и рассылки ведущей строки каждый процессор должен выполнить вычитание ведущей строки из каждой строки оставшейся части строк своей полосы матрицы A. Количество элементов строки, подлежащих обработке, также сокращается при исключении неизвестных, и текущее число элементов строки для вычислений оценивается величиной (n-i). Тем самым, сложность процедуры вычитания строк оценивается как 2x(n-i) операций (перед вычитанием ведущая строка умножается на масштабирующую величину aik/aii ). С учетом выполняемого количества итераций общее число операций параллельного варианта прямого хода метода Гаусса определяется выражением:$$T_p^1=\sum_{i=0}^{n-2} \left[ \frac{(n-i)}{p}+\frac{(n-i)}{p}\cdot 2(n-i) \right] = \frac1p \sum_{i=0}^{n-2} \left[ (n-i)+2(n-i)^2 \right] .$$

На каждой итерации обратного хода алгоритма Гаусса после рассылки вычисленного значения очередной неизвестной каждый процессор должен обновить значения правых частей для всех строк, расположенных на этом процессоре. Отсюда следует, что трудоемкость параллельного варианта обратного хода алгоритма Гаусса оценивается как величина:$$T_p^2 = \sum_{i=0}^{n-2} 2 \cdot (n-i)/p.$$

Просуммировав полученные выражения, можно получить$$T_p=\frac1p \sum_{i=0}^{n-2} [(n-i)+2(n-i)^2] + \frac2p \sum_{i=0}^{n-2} (n-i) = \frac1p \sum_{i=0}^{n-2} [3(n-i)+2(n-i)^2] = \frac1p \sum_{i=2}^{n} (3i+2i^2) .$$

Как результат выполненного анализа, показатели ускорения и эффективности параллельного варианта метода Гаусса могут быть определены при помощи соотношений следующего вида:$$S_p = \frac{(2n^3/3+n^2)}{\frac1p \sum\limits_{i=2}^{n}(3i+2i^2)}, \quad E_p = \frac{(2n^3/3+n^2)}{\sum\limits_{i=2}^{n}(3i+2i^2)} .$$

Полученные соотношения имеют достаточно сложный вид для оценивания. Вместе с тем можно показать, что сложность параллельного алгоритма имеет порядок ~(2n3/3)/p, и, тем самым, балансировка вычислительной нагрузки между процессорами в целом является достаточно равномерной.

Дополним сформированные показатели вычислительной сложности метода Гаусса оценкой затрат на выполнение операций передачи данных между процессорами. При выполнении прямого хода на каждой итерации для определения ведущей строки процессоры обмениваются локально найденными максимальными значениями в столбце с исключаемой переменной. Выполнение данного действия одновременно с определением среди собираемых величин наибольшего значения может быть обеспечено при помощи операции обобщенной редукции (функция MPI_Allreduce библиотеки MPI ). Всего для выполнения такой операции требуется log2p шагов, что с учетом количества итераций позволяет оценить время, необходимое для проведения операций редукции, при помощи следующего выражения:$$T_p^1(comm)=(n-1)\cdot\log_2p\cdot(\alpha+w/ \beta) ,$$ где, как и ранее, $$\alpha$$ – латентность сети передачи данных, $$\beta$$ – пропускная способность сети, w – размер пересылаемого элемента данных.

Далее также на каждой итерации прямого хода метода Гаусса выполняется рассылка выбранной ведущей строки. Сложность данной операции передачи данных:$$T_p^2(comm)=(n-1)\cdot\log_2p\cdot(\alpha+wn/ \beta) ,$$

При выполнении обратного хода алгоритма Гаусса на каждой итерации осуществляется рассылка между всеми процессорами вычисленного значения очередной неизвестной. Общее время, необходимое для выполнения подобных действий, можно оценить как:$$T_p^3(comm)=(n-1)\cdot\log_2p\cdot(\alpha+w/ \beta) ,$$

В итоге, с учетом всех полученных выражений, трудоемкость параллельного варианта метода Гаусса составляет:$$T_p = \frac1p \sum_{i=2}^n (3i+2i^2)\tau + (n-1)\cdot\log_2 p\cdot(3\alpha+w(n+2)/ \beta),$$ где $$\tau$$ есть время выполнения базовой вычислительной операции.

8.2.6. Программная реализация

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

1. Главная функция программы. Реализует логику работы алгоритма, последовательно вызывает необходимые подпрограммы.

// Программа 8.1. — Алгоритм Гаусса решения систем линейных уравнений
int ProcNum;       	// Число доступных процессоров 
int ProcRank;       	// Ранг текущего процессора
int *pParallelPivotPos;	// Номера строк, которые были выбраны ведущими
int *pProcPivotIter; 	// Номера итераций, на которых строки 
                   	// процессора использовались в качестве ведущих
void main(int argc, char* argv[]) {
  double* pMatrix;  		// Матрица линейной системы
  double* pVector;  		// Вектор правых частей линейной системы
  double* pResult;  		// Вектор неизвестных
  double *pProcRows; 	// Строки матрицы A
  double *pProcVector; 	// Блок вектора b
  double *pProcResult;	// Блок вектора x
  int     Size;       	// Размер матрицы и векторов
  int     RowNum;     	// Количество строк матрицы 
  double start, finish, duration;

  setvbuf(stdout, 0, _IONBF, 0);

  MPI_Init ( argc, argv );
  MPI_Comm_rank ( MPI_COMM_WORLD, ProcRank );
  MPI_Comm_size ( MPI_COMM_WORLD, ProcNum );
  
  if (ProcRank == 0)
    printf("Параллельный метод Гаусса для решения систем линейных уравнений\n");

  // Выделение памяти и инициализация данных
  ProcessInitialization(pMatrix, pVector, pResult, 
    pProcRows, pProcVector, pProcResult, Size, RowNum);

  // Распределение исходных данных
  DataDistribution(pMatrix, pProcRows, pVector, pProcVector, 
  Size, RowNum);

  // Выполнение параллельного алгоритма Гаусса
  ParallelResultCalculation(pProcRows, pProcVector, pProcResult, Size,
    RowNum);

  // Сбор найденного вектора неизвестных на ведущем процессе
  ResultCollection(pProcResult, pResult);
  
  // Завершение процесса вычислений
  ProcessTermination(pMatrix, pVector, pResult, pProcRows,
    pProcVector, pProcResult);

  MPI_Finalize();
}

Следует пояснить использование дополнительных массивов. Элементы массива pParallelPivotPos определяют номера строк матрицы, выбираемых в качестве ведущих, по итерациям прямого хода метода Гаусса. Именно в этом порядке должны выполняться затем итерации обратного хода для определения значений неизвестных системы линейных уравнений. Массив pParallelPivotPos является глобальным, и любое его изменение в одном из процессов требует выполнения операции рассылки измененных данных всем остальным процессам программы.

Элементы массива pProcPivotIter определяют номера итераций прямого хода метода Гаусса, на которых строки процесса использовались в качестве ведущих (т.е. строка i процесса выбиралась ведущей на итерации pProcPivotIter[i] ). Начальное значение элементов массива устанавливается нулевым, и, тем самым, нулевое значение элемента массива pProcPivotIter[i] является признаком того, что строка i процесса все еще подлежит обработке. Кроме того, важно отметить, что запоминаемые в элементах массива pProcPivotIter номера итераций дополнительно означают и номера неизвестных, для определения которых будут использованы соответствующие строки уравнения. Массив pProcPivotIter является локальным для каждого процесса.

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

Функция DataDistribution реализует распределение матрицы линейной системы и вектора правых частей между процессорами вычислительной системы.

Функция ResultCollection осуществляет сбор со всех процессов отдельных частей вектора неизвестных.

Функция ProcessTermination выполняет необходимый вывод результатов решения задачи и освобождает всю ранее выделенную память для хранения данных.

Реализация всех перечисленных функций может быть выполнена по аналогии с ранее рассмотренными примерами и предоставляется читателю в качестве самостоятельного упражнения.

2. Функция ParallelResultCalculation. Реализует логику работы параллельного алгоритма Гаусса, последовательно вызывает функции, выполняющие прямой и обратный ход метода.

// Функция для параллельного выполнения метода Гаусса
void ParallelResultCalculation (double* pProcRows, 
double* pProcVector, double* pProcResult, int Size, int RowNum) {
  ParallelGaussianElimination (pProcRows, pProcVector, Size,
  RowNum);
  ParallelBackSubstitution (pProcRows, pProcVector, pProcResult,
  Size, RowNum);
}

3. Функция ParallelGaussianElimination. Функция выполняет параллельный вариант прямого хода алгоритма Гаусса.

// Функция для параллельного выполнения прямого хода метода Гаусса
void ParallelGaussianElimination (double* pProcRows,
  double* pProcVector, int Size, int RowNum) {
  double MaxValue; 	// Значение ведущего элемента на процессоре
  int    PivotPos; 	// Положение ведущей строки в полосе линейной 
                  	// системы даннного процессора 
  // Структура для выбора ведущего элемента
  struct { double MaxValue; int ProcRank; } ProcPivot, Pivot;

  // pPivotRow используется для хранения ведущей строки матрицы и
  // соответствующего элемента вектора b
  double* pPivotRow = new double [Size+1];
  // Итерации прямого хода метода Гаусса
  for (int i=0; i<Size; i++) {
  
    // Нахождение ведущей строки среди строк процесса
    double MaxValue = 0;
    for (int j=0; j<RowNum; j++) {
      if ((pProcPivotIter[j] == -1)  
          (MaxValue < fabs(pProcRows[j*Size+i]))) {
        MaxValue = fabs(pProcRows[j*Size+i]);
        PivotPos = j;
      }
    }
    ProcPivot.MaxValue = MaxValue;
    ProcPivot.ProcRank = ProcRank;

    // Нахождение ведущего процесса (процесса, который содержит
    // максимальное значение переменной MaxValue)
    MPI_Allreduce(ProcPivot, Pivot, 1, MPI_DOUBLE_INT, 
MPI_MAXLOC, MPI_COMM_WORLD);

    // Рассылка ведущей строки
    if ( ProcRank == Pivot.ProcRank ){
      pProcPivotIter[PivotPos]= i; // номер итерации
      pParallelPivotPos[i]= pProcInd[ProcRank] + PivotPos;
    }
    MPI_Bcast(pParallelPivotPos[i], 1, MPI_INT, Pivot.ProcRank,
      MPI_COMM_WORLD); 

    if ( ProcRank == Pivot.ProcRank ){
      // заполнение ведущей строки
      for (int j=0; j<Size; j++) {
        pPivotRow[j] = pProcRows[PivotPos*Size + j];
      }
      pPivotRow[Size] = pProcVector[PivotPos];
    }
    MPI_Bcast(pPivotRow, Size+1, MPI_DOUBLE, Pivot.ProcRank, 
      MPI_COMM_WORLD);
  ParallelEliminateColumns(pProcRows, pProcVector, pPivotRow,
  Size, RowNum, i);
  }
}

Функция ParallelEliminateColumns проводит вычитание ведущей строки из строк процесса, которые еще не использовались в качестве ведущих (т.е. для которых элементы массива pProcPivotIter равны нулю).

4. Функция ParallelBackSubstitution. Функция реализует параллельный вариант обратного хода Гаусса.

// Функция для параллельного выполнения обратного хода метода Гаусса
void ParallelBackSubstitution (double* pProcRows, double* pProcVector,
  double* pProcResult, int Size, int RowNum) {
  int IterProcRank;	// Ранг процессора, который содержит ведущую строку
  int IterPivotPos;	// Положение ведущей строки в полосе процессора
  double IterResult;	// Вычисленное значение очередной неизвестной
  double val;

  // Итерации обратного хода метода Гаусса
  for (int i=Size-1; i>=0; i--) {

    // Вычисление ранга процессора, который содержит ведущую строку
    FindBackPivotRow(pParallelPivotPos[i], Size, IterProcRank,
      IterPivotPos);
    
    // Вычисление значения неизвестной
    if (ProcRank == IterProcRank) {
      IterResult = pProcVector[IterPivotPos]/pProcRows[IterPivotPos*Size+i];
      pProcResult[IterPivotPos] = IterResult;
    }
    // Рассылка значения очередной неизвестной
    MPI_Bcast(IterResult, 1, MPI_DOUBLE, IterProcRank, MPI_COMM_WORLD);

    // Обновление вектора правых частей
    for (int j=0; j<RowNum; j++) 
      if ( pProcPivotIter[j] < i ) {
        val = pProcRows[j*Size + i] * IterResult;
        pProcVector[j]=pProcVector[j] - val;
    }
  }
}

Функция FindBackPivotRow определяет строку, из которой можно вычислить значение очередного элемента результирующего вектора. Номер этой строки хранится в массиве pParallelPivotIter. По номеру функция FindBackPivotRow определяет номер процесса, на котором эта строка хранится, и номер этой строки в полосе pProcRows этого процесса.

8.2.7. Результаты вычислительных экспериментов

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

Эксперименты производились на вычислительном кластере Нижегородского университета на базе процессоров Intel Xeon 4 EM64T, 3000 МГц и сети Gigabit Ethernet под управлением операционной системы Microsoft Windows Server 2003 Standard x64 Edition и системы управления кластером Microsoft Compute Cluster Server (см. п. 1.2.3).

Для оценки длительности $$\tau$$ базовой скалярной операции проводилось решение системы линейных уравнений при помощи последовательного алгоритма и полученное таким образом время вычислений делилось на общее количество выполненных операций – в результате подобных экспериментов для величины $$\tau$$ было получено значение 4,7 нсек. Эксперименты, выполненные для определения параметров сети передачи данных, показали значения латентности a и пропускной способности b соответственно 47 мкс и 53,29 Мбайт/с. Все вычисления производились над числовыми значениями типа double, т.е. величина w равна 8 байт.

Результаты вычислительных экспериментов приведены в таблице 8.1. Эксперименты выполнялись с использованием двух, четырех и восьми процессоров.

Сравнение времени выполнения эксперимента $$T^*_p$$ и теоретической оценки Tp из (8.5) приведено в таблице 8.2 и на рис. 8.4.

8.3. Метод сопряженных градиентов

Рассмотрим теперь совершенно иной подход к решению систем линейных уравнений, при котором к искомому точному решению x* системы Ax=b строится последовательность приближенных решений x0, x1, ..., xk, ... При этом процесс вычислений организуется таким способом, что каждое очередное приближение дает оценку точного решения со все уменьшающейся погрешностью, и при продолжении расчетов оценка точного решения может быть получена с любой требуемой точностью. Подобные итерационные методы решения систем линейных уравнений широко используются в практике вычислений. К преимуществам итерационных методов можно отнести меньший объем (по сравнению, например, с методом Гаусса) необходимых вычислений для решения разреженных систем линейных уравнений, возможность более быстрого получения начальных приближений искомого решения, наличие эффективных способов организации параллельных вычислений. Дополнительная информация с описанием таких методов, рассмотрением вопросов сходимости и точности получаемых решений может быть получена, например, в [, ].

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

Напомним, что матрица А является симметричной, если она совпадает со своей транспонированной матрицей, т.е. А=АТ. Матрица А называется положительно определенной, если для любого вектора x справедливо: xTAx>0.

Известно (см., например, [, ]), что после выполнения n итераций алгоритма ( n есть порядок решаемой системы линейных уравнений ), очередное приближение xn совпадает с точным решением.

8.3.1. Последовательный алгоритм

Если матрица A симметричная и положительно определенная, то функция$$q(x)=\frac12 x^T \cdot A \cdot x - x^T b+c$$ имеет единственный минимум, который достигается в точке x*, совпадающей с решением системы линейных уравнений (8.2). Метод сопряженных градиентов является одним из широкого класса итерационных алгоритмов, которые позволяют найти решение (8.2) путем минимизации функции q(x).

Итерация метода сопряженных градиентов состоит в вычислении очередного приближения к точному решению в соответствии с правилом:$$x^k = x^{k-1} + s^k d^k$$

Тем самым, новое значение приближения xk вычисляется с учетом приближения, построенного на предыдущем шаге xk-1, скалярного шага sk и вектора направления dk.

Перед выполнением первой итерации векторы x0 и d0 полагаются равными нулю, а для вектора g0 устанавливается значение, равное - b. Далее каждая итерация для вычисления очередного значения приближения xk включает выполнение четырех шагов:

Шаг 1: Вычисление градиента:$$g^k=A\cdot x^{k-1}-b;$$

Шаг 2: Вычисление вектора направления:$$d^k=-g^k+\frac{((g^k)^T,g^k)}{((g^{k-1})^T,g^{k-1})}d^(k-1),$$ где (gT,g) есть скалярное произведение векторов.

Шаг 3: Вычисление величины смещения по выбранному направлению:$$s^k=\frac{(d^k,g^k)}{(d^k)^T\cdot A\cdot d^k}$$

Шаг 4: Вычисление нового приближения:$$x^k = x^{k-1}+s^k d^k$$

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

t1=2n2+13n.

Как уже отмечалось ранее, для нахождения точного решения системы линейных уравнений с положительно определенной симметричной матрицей необходимо выполнить n итераций. Таким образом, для нахождения решения системы необходимо выполнить$$T_1=2n^3+13n^2,$$ и, тем самым, сложность алгоритма имеет порядок O(n3).

Поясним выполнение метода сопряженных градиентов на примере решения системы линейных уравнений вида:

3x0 -x1  = 3,
-x0 +3x1 = 7.

Эта система уравнений второго порядка обладает симметричной положительно определенной матрицей, для нахождения точного решения этой системы достаточно провести всего две итерации метода.

На первой итерации было получено значение градиента g1=(-3,-7), значение вектора направления d1=(3, 7), значение величины смещения s1=0,439. Соответственно, очередное приближение к точному решению системы x1=(1,318, 3,076).

На второй итерации было получено значение градиента g2=(-2,121, 0,909), значение вектора направления d2=(2,397, -0,266), а величина смещения – s2=0,284. Очередное приближение совпадает с точным решением системы x2=(2, 3).

На рис. 8.5 представлена последовательность приближений к точному решению, построенная методом сопряженных градиентов (в качестве начального приближения x0 выбрана точка (0,0) ).

8.3.2. Организация параллельных вычислений

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

Анализ соотношений (8.8) – (8.11) показывает, что основные вычисления, выполняемые в соответствии с методом, состоят в умножении матрицы A на векторы x и d, и, как результат, при организации параллельных вычислений может быть полностью использован материал, изложенный в лекции 6.

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

8.3.3. Анализ эффективности

Выберем для дальнейшего анализа эффективности получаемых параллельных вычислений параллельный алгоритм матрично-векторного умножения при ленточном горизонтальном разделении матрицы и при полном дублировании всех обрабатываемых векторов.

Трудоемкость последовательного метода сопряженных градиентов была уже определена ранее в (8.12).

Определим время выполнения параллельной реализации метода сопряженных градиентов. Вычислительная сложность параллельных операций умножения матрицы на вектор при использовании схемы ленточного горизонтального разделения матрицы составляет:$$T_p^1(calc)=2n\lceil n/p\rceil\cdot(2n-1)$$ (см. лекцию 6).

Как результат, при условии дублирования всех вычислений над векторами общая вычислительная сложность параллельного варианта метода сопряженных градиентов равна:$$T_p(calc)=n(2\cdot\lceil n/p\rceil\cdot(2n-1)+13n).$$

С учетом полученных оценок показатели ускорения и эффективности параллельного алгоритма могут быть выражены при помощи соотношений:$$S_p=\frac{2n^3+13n^2}{n(2\lceil n/p\rceil\cdot(2n-1)+13n)}, \quad E_p=\frac{2n^3+13n^2}{p\cdot n(2\lceil n/p\rceil\cdot(2n-1)+13n)}$$

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

Уточним теперь приведенные выражения – учтем длительность выполняемых вычислительных операций и оценим трудоемкость операции передачи данных между процессорами. Как можно заметить, информационное взаимодействие процессоров возникает только при выполнении операций умножения матрицы на вектор. С учетом результатов лекции 6 коммуникационная сложность рассматриваемых параллельных вычислений равна:$$T_p(comm)=2n(\alpha\cdot\lceil\log p\rceil+w(n/p)(p-1)/ \beta),$$ где $$\alpha$$ – латентность, $$\beta$$ – пропускная способность сети передачи данных, а w есть размер элемента упорядочиваемых данных в байтах.

Окончательно, время выполнения параллельного варианта метода сопряженных градиентов для решения систем линейных уравнений принимает вид:$$T_p=n\cdot[(2\lceil n/p\rceil\cdot(2n-1)+13n)\cdot\tau+2(\alpha\cdot\lceil\log p\rceil+w(n/p)(p-1)/ \beta)]$$ где $$\tau$$ есть время выполнения базовой вычислительной операции.

8.3.4. Результаты вычислительных экспериментов

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

Результаты вычислительных экспериментов приведены в таблице 8.3. Эксперименты проводились на вычислительных системах, состоящих из двух, четырех и восьми процессоров.

Сравнение времени выполнения эксперимента $$T^*_p$$ и теоретической оценки Tp из (8.13) приведено в таблице 8.4 и на рис. 8.7.

8.4. Краткий обзор лекции

Данная лекция посвящена проблеме параллельных вычислений при решении систем линейных уравнений. Изложение учебного материала проводится с использованием двух широко известных алгоритмов: метода Гаусса как примера прямого алгоритма решения задачи и итерационного метода сопряженных градиентов.

Параллельный вариант метода Гаусса (подраздел 8.2) основывается на ленточном разделении матрицы между процессорами с использованием циклической схемы распределения строк, что позволяет сбалансировать вычислительную нагрузку. Для определения параллельного варианта метода проводится полный цикл проектирования – определяются базовые подзадачи, выделяются информационные взаимодействия, обсуждаются вопросы масштабирования, выводятся оценки показателей эффективности, предлагается схема программной реализации и приводятся результаты вычислительных экспериментов. Как видно из графика, представленного на рис. 8.8, параллельный алгоритм Гаусса демонстрирует достаточно высокие показатели ускорения и эффективности.

Важный момент при рассмотрении параллельного варианта метода сопряженных градиентов (подраздел 8.3) состоит в том, что способ параллельных вычислений для этого метода может быть получен через параллельные алгоритмы выполняемых вычислительных действий – операций умножения матрицы на вектор, скалярного произведения векторов, сложения и вычитания векторов. Выбранный для учебного примера подход состоит в распараллеливании вычислительно наиболее трудоемкой операции умножения матрицы на вектор, в то время как все вычисления над векторами дублируются на каждом процессоре для уменьшения числа выполняемых операций передачи данных — показатели эффективности подобного способа организации параллельных вычислений представлены на рис. 8.8.

8.5. Обзор литературы

Проблема численного решения систем линейных уравнений широко рассматривается в литературе. Для учебного рассмотрения могут быть рекомендованы работы [, , , ]. Широкое обсуждение вопросов параллельных вычислений для решения данного типа задач выполнено в работах [, ].

При рассмотрении вопросов программной реализации параллельных методов может быть рекомендована работа []. В ней рассматривается хорошо известная и широко используемая в практике параллельных вычислений программная библиотека численных методов ScaLAPACK.

8.6. Контрольные вопросы

8.7. Задачи и упражнения