5.2 Epetra: представление векторов и матриц

Вектора и матрицы являются фундаментальными структурами данных, которые используются почти во всех численных методах. В рамках библиотеки Trilinos при создании векторов, как правило, используются классы из пакета Epetra.

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

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

Рассмотрим классы, реализующие сущность коммуникатора. Виртуальный класс Epetra_Comm предоставляет интерфейс, инкапсулирующий общую информацию и методы, необходимые для организации и запуска других классов на компьютерах с параллельной или последовательной архитектурой. Создание объекта класса Epetra_Comm является необходимым при создании всех объектов класса Epetra_Map, который в свою очередь необходим для всех других классов из Epetra.

Epetra_Comm имеет две основные реализаций:

Большинство методов класса Epetra_Comm подобны функциям MPI. В классе имеются такие методы как MyPID(), NumProc(), Barrier(), Broadcast(), SumAll(), GatherAll(), MaxAll(), MinAll(), ScanSum(). Например, число процессов в коммуникаторе NumProc и идентификатор текущего процесса MyPID могут быть получены следующим образом:

int NumProc = Comm.NumProc(); 
int MyPID = Comm.MyPID();

Класс Epetra_Map является множеством целочисленных элементов, распределенным между процессами и содержащим следующую информацию:

Существуют три способа создания карты (map). Проще всего указать общее (глобальное) число элементов, а Epetra сам решит, как распределять между процессами:

Epetra_Map Map(NumGlobalElements, 0, Comm);

В этом случае конструктор класса принимает глобальную размерность вектора и объект Epetra_Comm. В результате, для каждого процесса будет определено свое локальное множество элементов.

При втором способе создания конструктору объекта передаем локальное количество элементов:

Epetra_Map Map(-1, NumMyElements, 0, Comm);

В этом случае создатся вектор размера i=0NumProc-1NumMyElements.

Третий способ создания Epetra_Map заключается в определении на каждом процессе двух параметров: количества локальных элементов и глобальных индексов для каждого локального элемента. Рассмотрим данный способ на примере, в котором создается вектор из 5 элементов, распределенный между двумя процессами p0 и p1. Процесс р0 владеет элементами 0 и 4, а процесс p1 – элементами 1, 2 и 3.

#include "Epetra_Map.h" 
// ... 
MyPID = Comm.MyPID(); 
switch( MyPID ) { 
case 0: 
    MyElements = 2; 
    MyGlobalElements = new int[MyElements]; 
    MyGlobalElements[0] = 0; 
    MyGlobalElements[1] = 4; 
    break; 
case 1: 
    MyElements = 3; 
    MyGlobalElements = new int[MyElements]; 
    MyGlobalElements[0] = 1; 
    MyGlobalElements[1] = 2; 
    MyGlobalElements[2] = 3; 
    break; 
} 
Epetra_Map Map(-1, MyElements, MyGlobalElements, 0, Comm);

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

int NumGlobalElements = Map.NumGlobalElements(); 
int NumMyElements = Map.NumMyElements();

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

int * MyGlobalElements = Map.MyGlobalElements();

которая возвращает указатель на вектор глобальной индексации.

Следует заметить, что классы в Epetra перегружают оператор , что делает простым вывод информации про класс.

Для создания последовательного вектора, содержащего double элементы, используется класс SerialDenseVector:

#include "Epetra_SerialDenseVector.h" 
Epetra_SerialDenseVector DoubleVector(Length);

Подобным образом можно создать последовательный вектор из integer значений

#include "Epetra_IntSerialDenseVector.h" 
Epetra_SerialIntDenseVector IntVector(Length);

Рекомендуется создавать вектора с использованием указателей в динамической памяти (используя new), поскольку при создании вектора в стеке он автоматически удалится при выходе из области видимости.

Для обращения к элементам вектора используются операторы [] и (). Оба оператора возвращают ссылки на элемент вектора.

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

Epetra_Vector x(Map);

При вызове конструктора Epetra_Vector выделяется память и задаются нулевые значения элементов.

Создание вектора может происходить посредством вызова конструктора копий:

Epetra_Vector y(x);

Для задания значения конкретному элементу используем оператор [], например:

x[i] = 1.0*i;

где значение задается i-му элементу вектора x.

Матрицы в Epetra, как и вектора, могут быть как последовательными, так и распределенными. При создании последовательной матрицы используется класс Epetra_SerialDenseMatrix. Создание плотной матрицы D размерности n на m:

#include "Epetra_SerialDenseMatrix.h" 
Epetra_SerialDenseMatrix D(n,m);

Также можно создать матрицу нулевой размерности:

Epetra_SerialDenseMatrix D();

Затем можно задать размерность уже созданной матрице:

D.Shape(n,m);

Для получения доступа к элементу расположенному в i-ой строке j-го столбца используется оператор ()

A(i,j)

или оператор []

A[j][i]

Приведем пример использования рассмотренного класса:

int NumRowsA = 2, NumColsA = 2; 
int NumRowsB = 2, NumColsB = 1; 
Epetra_SerialDenseMatrix A, B; 
A.Shape(NumRowsA, NumColsA); 
B.Shape(NumRowsB, NumColsB); 
Epetra_SerialDenseMatrix AtimesB; 
AtimesB.Shape(NumRowsA,NumColsB); 
double alpha = 1.0, beta = 1.0; 
AtimesB.Multiply(N,N,alpha, A, B, beta); 
cout << AtimesB;

В примере описан код умножения квадратных матриц A и B.

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

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

В качестве примера рассмотрим создание распределенной разреженной матрицы, получаемой в результате конечно-разностной аппроксимации оператора Лапласа

⌊ 2    - 1    0     ...    0⌋
||- 1    2    - 1    ...    0||
| ..     ..     ..     ..     ..|
⌈ .     .     .     .     .⌉
  0     0     0     ...    2

В начале определим глобальную размерность

int NumGlobalElements = 5;

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

Epetra_Map Map(NumGlobalElements,0,Comm); 
int NumMyElements = Map.NumMyElements(); 
int * MyGlobalElements = Map.MyGlobalElements( );

Потом можем создать матрицу

Epetra_CrsMatrix A(Copy, Map, 3);

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

Или матрицу можно создать по-другому, учитывая, что в первой и последней строке матрицы содержится только два ненулевых элемента

int * NumNz = new int[NumMyElements]; 
for( int i=0 ; i<NumMyElements ; i++ ) { 
    if( MyGlobalElements[i]==0 ||   MyGlobalElements[i] == NumGlobalElements-1) 
        NumNz[i] = 2; 
    else 
        NumNz[i] = 3; 
} 
Epetra_CrsMatrix A(Copy, Map, NumNz);

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

double * Values = new double[2]; 
Values[0] = -1.0; 
Values[1] = -1.0; 
int * Indices = new int[2]; 
double two = 2.0; 
int NumEntries; 
for( int i=0 ; i<NumMyElements; ++i ) { 
    if (MyGlobalElements[i]==0) { 
        Indices[0] = 1; 
        NumEntries = 1; 
    } else if (MyGlobalElements[i] == NumGlobalElements-1) { 
        Indices[0] = NumGlobalElements-2; 
        NumEntries = 1; 
    } else { 
        Indices[0] = MyGlobalElements[i]-1; 
        Indices[1] = MyGlobalElements[i]+1; 
        NumEntries = 2; 
    } 
    A.InsertGlobalValues(MyGlobalElements[i], NumEntries,Values, Indices); 
    A.InsertGlobalValues(MyGlobalElements[i], 1, &two, MyGlobalElements+i); 
}

И в конце мы вызываем следующую функцию:

A.FillComplete();

В пакете Epetra имеется также класс для работы со временем Epetra_Time. Для измерения времени выполнения какого-либо куска кода в вашей программе нужно сначала получить текущее время

Epetra_Time time(Comm);

а затем вызвать

time.ResetStartTime();