6.2 Представление векторов и матриц

Вектор – это простейший объект в PETSC. Вектора используются для хранения решения, правой части систем линейных уравнений и т.п. [11], [3]

В PETSC есть два основных типа вектора: последовательный и параллельный (основанный на MPI). Чтобы создать последовательный вектор с m компонентами, вызываем следующую функцию:

    VecCreateSeq(PETSC_COMM_SELF, int m, Vec *x);

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

    VecCreateMPI(MPI_Comm comm, int m, int M, Vec *x);

создает вектор, который распределен по всем процессам в коммуникаторе comm, где m указывает количество компонентов, которые нужно хранить для локального процесса, а M есть общая размерность вектора. Локальную или глобальную размерность вектора можно задать как PETSC_DECIDE, тогда PETSc сам определит одну из размерностей. В общем случае можно использовать функцию

    VecCreate(MPI_Comm comm, int m, int M, Vec *v); 
    VecSetFromOptions(v);

которая автоматически генерируют соответствующий векторный тип (последовательный или параллельный) для всех процессов в коммуникаторе comm.

Задать одно и тоже значение для всех компонент вектора можно с помощью функции

    VecSet(Scalar *value, Vec x);

Чтобы компонентам вектора присвоить различные значения вызываем

    VecSetValues(Vec x, int n, int *indices, Scalar *values, INSERT_VALUES);

где n – количество присваиваемых компонент, массив indices содержит индексы глобальных компонент, а значения должны лежать в массиве values.

Когда все значения уже размещены с помощью VecSetValues(), можно вызвать

    VecAssemblyBegin(Vec x); 
    VecAssemblyEnd(Vec x);

Вектор также можно создать с использованием дублирования существующего вектора

    VecDuplicate(Vec old, Vec *new);

Для того чтобы вывести вектор на экран вызываем функцию

    VecView(Vec x, PetscViewer v);

где в качестве вьювера задаем VIEWER_STDOUT_WORLD.

Когда вектор больше не используется его следует удалить вызвав функцию:

    VecDestroy(Vec x);

В таблице 6.1 представлены основные операции для работы с векторами.




Название функции Операция


VecAXPY(Scalar *a, Vec x, Vec y); y = y + a * x


VecAYPX(Scalar *a,Vec x, Vec y); y = x + a * y


VecWAXPY(Scalar *a,Vec x,Vec y, Vec w); w = a * x + y


VecAXPBY(Scalar *a,Scalar *,Vec x,Vec y); y = a * x + b * y


VecScale(Scalar *a, Vec x); x = a * x


VecDot(Vec x, Vec y, Scalar *r); r = x′* y


VecTDot(Vec x, Vec y, Scalar *r); r = x′* y


VecNorm(Vec x,NormType type, double *r); r = ||x||type


VecSum(Vec x, Scalar *r); r = xi


VecCopy(Vec x, Vec y); y = x


VecSwap(Vec x, Vec y); y = x while x = y


VecPointwiseMult(Vec x, Vec y, Vec w); wi = xi * yi


VecPointwiseDivide(Vec x, Vec y, Vec w); wi = xi∕yi


VecMDot(int n,Vec x, Vec *y,Scalar *r); r[i] = x′* y[i]


VecMTDot(int n,Vec x, Vec *y,Scalar *r); r[i] = x′* y[i]


VecMAXPY(int n, Scalar *a,Vec y, Vec *x);y = y + ai * x[i]


VecMax(Vec x, int *idx, double *r); r = maxxi


VecMin(Vec x, int *idx, double *r); r = minxi


VecAbs(Vec x); xi = |xi|


VecReciprocal(Vec x); xi = 1∕xi


VecShift(Scalar *s,Vec x); xi = s + xi



Таблица 6.1: Основные векторные операции в PETSc

Использование матриц в PETSc похоже на использование векторов. Пользователь может создать новую параллельную или последовательную матрицу A, у которой M строк и N столбцов, также может указать каждому процессу локальное число строк и столбцов через параметры m и n. Для создания вызывается

    MatCreate (MPI_Comm comm,int m,int n,int M,int N, Mat* A);

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

Значения могут быть установлены командой:

    MatSetValues (Mat A,int m,int *im,int n,int *in, 
                  PetscScalar *values, INSERT_VALUES);

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

    MatAssemblyBegin (Mat A,MAT_FINAL_ASSEMBLY); 
    MatAssemblyEnd (Mat A,MAT_FINAL_ASSEMBLY);

В таблице 6.2 представлены основные матричные операции в PETSc.




Название функции Операция


MatAXPY(Mat Y, PetscScalar a, Mat X, MatStructure str); Y = Y + a * X


MatMult(Mat A,Vec x, Vec y); y = A * x


MatMultAdd(Mat A,Vec x, Vec y,Vec z); z = y + A * x


MatMultTranspose(Mat A,Vec x, Vec y); y = AT * x


MatMultTransposeAdd(Mat A,Vec x, Vec y, Vec z); z = y + AT * x


MatNorm(Mat A,NormType type, double *r); r = ||A||type


MatDiagonalScale(Mat A,Vec l,Vec r); A = diag(l) * A * diag(r)


MatScale(Mat A,PetscScalar a); A = a * A


MatConvert(Mat A,MatType type,Mat *B); B = A


MatCopy(Mat A,Mat B,MatStructure); B = A


MatGetDiagonal(Mat A,Vec x); x = diag(A)


MatTranspose(Mat A,MatReuse,Mat* B); B = AT


MatZeroEntries(Mat A); A = 0


MatShift(Mat Y,PetscScalar a); Y = Y + a * I



Таблица 6.2: Основные операции с матрицами в PETSc

По умолчанию матрицы в PETSc представлены в общем разреженном AIJ формате (так называемый the Yale sparse matrix format или compressed sparse row format, CSR). Рассмотрим некоторые детали эффективного использования этого формата для приложений большого размера.

Последовательные AIJ разреженные матрицы. В формате AIJ хранятся ненулевые элементы по строкам вместе с сопутствующим массивом соответствующих номеров столбцов и массивом указателей на начало каждой строки. Диагональные элементы матрицы хранятся с остальными ненулевыми элементами (не отдельно). Чтобы создать последовательную AIJ матрицу А с m строками и n столбцами вызываем функцию

    MatCreateSeqAIJ(PETSC_COMM_SELF, int m, int n, int nz, int *nnz, Mat *A);

где nz или nnz могут быть использованы для предварительного распределения матричной памяти. Пользователь может установить nz = 0 и nnz = PETSC_NULL, чтобы управлять всем размещением матричной памяти.

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

    int nnz[m]; 
    nnz[0] = <nonzeros in row 0> 
    nnz[1] = <nonzeros in row 1> 
    ... 
    nnz[m-1] = <nonzeros in row m-1>

Указав опцию -log_info при запуске вашей программы можно получить информацию о сборке матрицы. Рассмотрим следующий пример для матричного формата MATSEQ AIJ:

    MatAssemblyEnd_SeqAIJ:Matrix size 10 X 10; storage space:20 unneeded, 100 used 
    MatAssemblyEnd_SeqAIJ:Number of mallocs during MatSetValues is 0

Первая строка показывает, что пользователь указал 3000 областей памяти, но только 1000 областей была использована. Вторая строка отмечает, что пользователь указал достаточную область, чтобы PETSc не требовал дополнительной памяти для размещения. В следующем примере пользователь не указывает достаточно пространства.

    MatAssemblyEnd_SeqAIJ:Matrix size 10 X 10; storage space:47 unneeded, 1000 used 
    MatAssemblyEnd_SeqAIJ:Number of mallocs during MatSetValues is 40000

Параллельные AIJ разреженные матрицы. Параллельные разреженные матрицы с AIJ форматом могут быть созданы командой

    MatCreateMPIAIJ(MPI_Comm comm, int m, int n, int M, int N, 
                    int d_nz,int *d_nnz, int o_nz, int *o_nnz, Mat *A);

Аргумент А есть вновь созданная матрица, а аргументы m, n, M и N указывают число локальных и глобальных строк и столбцов соответственно.

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

Распределенные массивы (Distributed arrays – DAs) используются совместно с векторами PETSc и предназначены для работы с логическими регулярными прямоугольными сетками (структурированными), когда необходим обмен нелокальными данными перед некоторыми локальными вычислениями. Эти массивы спроектированы только для случая, когда данные можно представить хранящимися в многомерном стандартном массиве, следовательно, DA предназначен для обмена векторной информации и не предназначен для хранения матриц. Объект DA используется для управления параллельными коммуникациями, поскольку он содержит информацию о сетке и распределении ее по процессам.

Например, типичная ситуация, которая встречается при параллельном решении ДУЧП состоит в том, что для оценки локальной функции f(x) каждому процессору требуется его локальная часть вектора х вместе с теневыми точками (Граничные части вектора, принадлежащие соседним процессорам). На рис. 6.2 представлены теневые точки для шестого процессора в двумерной регулярной параллельной сетке. Каждый прямоугольник представляет процессор, теневые точки показаны серыми.


PIC


Рис. 6.2: Теневые точки


Для создания двухмерной структуры данных используется следующая функция:

    DACreate2d(MPI_Comm comm, DAPeriodicType wrap, 
               DAStencilType st, int M, int N,int m, int n, 
               int dof, int s, int *lx, int *ly, DA *da);

Здесь M и N – глобальное(общее) количество точек по каждому направлению, m и n – количество процессов по каждому направлению (m *n – общее количество запущенных процессов), s – тип коммуникации, который определяет шаблон аппроксимации (DA_STENCIL_STAR – пятиточечный шаблон или DA_STENCIL_BOX – девятиточечный шаблон).

Создание трехмерной структуры выглядит следующим образом:

    DACreate3d(MPI_Comm comm, DAPeriodicType wrap, 
               DAStencilType stencil_type, 
               int M, int N, int P, int m, int n, int p, 
               int w, int s, int *lx, int *ly, int *lz, DA *inra);

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

    DACreateGlobalVector(DA da, Vec *g); 
    DACreateLocalVector(DA da, Vec *l);

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

    DAGlobalToLocalBegin(DA da, Vec g, InsertMode iora, Vec l); 
    DAGlobalToLocalEnd(DA da, Vec g, InsertMode iora, Vec l);

Эти вектора будут служить строительными блоками для локального и глобального решения ДУЧП и др. Если необходимы дополнительные вектора с таким же распределением по процессам, они могут быть получены дублированием с использованием функций, VecDuplicate() или VecDuplicateVecs().

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

    DAGetLocalVector(DA da, Vec *l); 
    /* Use local vector l */ 
    ... 
    DARestoreLocalVector(DA da, Vec *l);

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

    DAVecGetArray(DA da, Vec l, void *array); 
    /* Work with array */ 
    ... 
    DAVecRestoreArray(DA da, Vec l, void *array);

где array – многомерный массив.

Глобальные индексы левого нижнего угла локальной части массива могут быть получены командами

    DAGetCorners(DA da, int *x, int *y, int *z, int *m, int *n, int *p);

В PETSc имеются модули для работы с неструктурированными сетками, но мы их здесь не рассматриваем.