Часть III.

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

Содержание

1.Введение.

Настоящее пособие представляет собой третью часть методических материалов, предназначенных для пользователей Комплекса программ для решения задач линейной алгебры с использованием распределенной памяти, реализуемого в НИВЦ МГУ (http://www.srcc.msu.su/num_anal/par_prog/).

Указанный Комплекс программ ориентирован на реализацию программ на ФОРТРАН'е - 77 в стиле SPMD (Single Program Multiple Data) с использованием для обмена информацией между параллельными процессами передачи сообщений с помощью примитивов MPI (Massage Passing Interface).

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

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

Повторим здесь лишь вкратце основные действия, требующиеся от пользователя программ Комплекса.
Прежде чем обратиться к какой - либо из параллельных программ пользователь должен выбрать для ее выполнения и описать так называемую "решетку процессов", т.е. "выстроить" параллельные процессы в общем случае в двумерное упорядоченное множество, состоящее из NPROW строк, в каждой из которых содержится по NPCOL процессов. Таким образом общее число используемых процессов  NP = NPROW * NPCOL.

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

По окончании вычислений необходимо также сделать выход из определенного ранее контекста (т.е. освободить использованные процессы решетки).

Кроме выбора решетки процессов, необходимо распределить исходные матрицы (векторы), обрабатываемые программами по параллельным процессам выбранной решетки.

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

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

В связи с вышесказанным, при вызове подпрограмм Комплекса помимо адреса локальной части матрицы (вектора) им необходимо передать в качестве фактического параметра адрес "дескрипторного массива" или "дескриптора". Информация, содержащаяся в дескрипторе, позволяет всегда однозначно определить в локальной памяти, какого процесса и на каком именно месте располагается каждый элемент глобальной матрицы  A ( I, J ) (см. п.5.2).

Действия по инициализации решетки процессов и определению контекста выполняются подпрограммами пакета BLACS.

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

В п.7 настоящей части пособия достаточно подробно рассматриваются вопросы обработки ошибок и выдачи диагностических сообщений.

2. Краткое описание содержания раздела Комплекса для решения
линейной симметричной проблемы собственных значений.

Пусть  A является вещественной симметричной или эрмитовой матрицей порядка  n.

Симметричная проблема собственных значений (SEP) состоит в нахождении собственных значений  l и соответствующих собственных векторов  z, не равных нулю, удовлетворяющих уравнению 

       A z  =  l z ,       A  =  AT .

Для эрмитовой проблемы собственных значений 

       A z  =  l z ,       A  =  AH .

В обоих случаях собственные значения  l являются вещественными.

Когда требуется вычислить все собственные значения и все собственные векторы, мы представляем матрицу  A в виде: 

        A  =  Z L ZT ,   где

L - диагональная матрица, диагональные элементы которой являются
      собственными значениями,

Z - ортогональная (или унитарная) матрица, столбцы которой являются
      собственными векторами.

Это - классическая спектральная факторизация матрицы  A.

Подробнее, вычисление состоит из следующих этапов.

2.1. Вещественная симметричная или эрмитова матрица  A приводится к вещественной трехдиагональной форме T.

Если  A является вещественной симметричной, разложение будет следующим: 

 
      A  =  Q T QT ,      где  Q - ортогональная матрица,
                                       а  T - симметричная трехдиагональная.

Если  A является эрмитовой матрицей, разложение имеет вид 

      A  =  Q T QH ,      где  Q - унитарная матрица,
                                        а  T - вещественная симметричная трехдиагональная.

2.2. Далее вычисляются собственные значения и собственные векторы вещественной симметричной трехдиагональной матрицы  T.

Если вычисляются все собственные значения и собственные векторы, этот процесс эквивалентен факторизации  T следующим образом: 

      T  =  S L ST ,      где  S - ортогональная матрица,
                                         L - диагональная матрица.

Диагональные элементы матрицы L являются собственными значениями матрицы  T, которые также являются собственными значениями матрицы  A, а столбцы матрицы  S являются собственными векторами матрицы  T; собственными векторами матрицы  A являются столбцы матрицы Z = Q S, т.к. A = Z L ZT (или A = Z L ZH, когда A - эрмитова).

В вещественном случае разложение матрицы  A выполняется базовой подпрограммой PDSYTRD. В комплексном случае - базовой подпрограммой PZHETRD. Подпрограммы PDSYTRD (или PZHETRD) представляют матрицу  Q как произведение элементарных матриц отражения.

Когда требуется вычислить все собственные значения и собственные векторы матрицы T (в целевых подпрограммах PDSYEV1, PZHEEV1, PDSYEV2, PZHEEV2) используется подпрограмма DSTEQR2 (ZSTEQR2). Эта подпрограмма является модификацией версии подпрограммы DSTEQR (ZSTEQR) из пакета LAPACK [12]. Она выполняет меньшее число итераций, чем подпрограмма пакета LAPACK, благодаря усовершенствованию техники предсказания (the look - ahead technique). Некоторые дополнительные модификации позволяют каждому процессу выполнять частичную корректировку матрицы  Q. Эта подпрограмма вычисляет собственные значения и собственные векторы симметричной трехдиагональной матрицы, используя неявные  QL или  QR методы.

Когда требуется вычислить некоторые или все собственные значения (в целевых подпрограммах PDSYEV3, PZHEEV3, PDSYEV4, PZHEEV4, PDSYEV5, PZHEEV5, PDSYEV6, PZHEEV6) используется базовая подпрограмма PDSTEBZ. Она позволяет вычислять некоторые собственные значения либо лежащие в вещественном интервале, либо собственные значения, принадлежащие интервалу индексов от  i - ого до  j - ого в массиве собственных значений, упорядоченных по возрастанию. Вычисление может производиться с высокой точностью. Либо можно сделать вычисление более быстрым, но с достижением более низкой точности.

Базовая подпрограмма PDORMTR (или PZUNMTR - в комплексном случае) обеспечивает умножение некоторой матрицы на матрицу  Q без ее (Q) явного формирования, что используется для обратной трансформации собственных векторов матрицы  T в собственные векторы матрицы  A. Это выполняется базовой подпрограммой PDSTEIN.

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

Без переортогонализации обратная итерация может выдавать векторы, которые дают большие скалярные произведения. Для исправления этого большинство реализаций обратной итерации в таких подпрограммах, как DSTEIN (из пакета LAPACK), выполняют переортогонализацию в тех случаях, когда собственные значения различаются меньше, чем на  10 - 3 || T ||. В результате собственные векторы, вычисленные подпрограммой DSTEIN, почти всегда являются ортогональными, но затраты возрастают до величины  O (n3)  операций.

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

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

PDSYEV1 -
PZHEEV1    		Вычисление всех собственных значений вещественной симметричной (эрмитовой) матрицы с использованием  QL или QR - алгоритма
PDSYEV2 -
PZHEEV2    		Вычисление всех собственных значений и собственных векторов вещественной симметричной (эрмитовой) матрицы с использованием  QL или QR - алгоритма
PDSYEV3 -
PZHEEV3    		Вычисление собственных значений вещественной симметричной (эрмитовой) матрицы, принадлежащих заданному интервалу с использованием  QL или QR - алгоритма
PDSYEV4 -
PZHEEV4    		Вычисление собственных значений, принадлежащих заданному интервалу, и соответствующих собственных векторов вещественной симметричной (эрмитовой) матрицы с использованием  QL или QR - алгоритма
PDSYEV5 -
PZHEEV5    		Вычисление собственных значений, принадлежащих заданному интервалу индексов, вещественной симметричной (эрмитовой) матрицы с использованием  QL или QR - алгоритма
PDSYEV6 -
PZHEEV6    		Вычисление собственных значений вещественной симметричной (эрмитовой) матрицы, принадлежащих заданному интервалу индексов, и соответствующих собственных векторов с использованием  QL или QR - алгоритма

В Приложении 3_1 приводится конкретный пример описания целевой подпрограммы для вычисления всех собственных значений и собственных векторов вещественной симметричной матрицы;
а в Приложении 3_2 - пример обращения к целевой подпрограмме из Приложения 3_1 для случая вещественной симметричной матрицы, который подробно рассматривается в п.5.

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

Кроме того, в Приложении 3_3 приводится пример составления головной программы при чтении исходной матрицы из файла и записи результатов в файл.

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

Пусть  A и  B являются вещественными симметричными или эрмитовыми матрицами порядка  n, а матрица  B является также положительно определенной.

Обобщенная симметричная проблема собственных значений (GSEP) состоит в нахождении собственных значений  l  и соответствующих собственных векторов  z,  не равных нулю, удовлетворяющих одному из видов уравнений. 

(1)           A z   =   l B z ,

(2)         A B z  =  l z ,

(3)         B A z  =  l z .

При этом для всех видов уравнений (1), (2) или (3) искомые значения  l  всегда являются вещественными.

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

    ZT A Z   =   L   (или  ZH A Z  =  L  -  для эрмитовых матриц),
                             при решении уравнений вида (1) или (3),

  Z-1 A Z-T  =  I     (или  Z-1 A Z-H = I  -  для эрмитовых матриц),
                              при решении уравнения вида (2),

где  L - диагональная матрица с расположенными на диагонали
             собственными значениями.

Матрицы  Z  удовлетворяют также уравнениям 

     ZT B Z  =  I   (или  ZH B Z  =  I  -  для эрмитовых матриц),
                          при решении уравнений вида (1) или (2),
  ZT B-1 Z  =  I  (или  ZH B-1 Z  =  I  -  для эрмитовых матриц),
                          при решении уравнения вида (3).

Подробнее, вычисление состоит из следующих этапов.

3.1. Каждый из указанных видов обобщенной проблемы собственных значений (1), (2) или (3) может быть приведен к стандартной (линейной) симметричной проблеме собственных значений, посредством использования разложения Холецкого для матрицы  B:

     B = L LT   или   B = UT U  -  для вещественного симметричного случая,
     B = L LH   или   B = UH U  -  для случая эрмитовых матриц.

Тогда при   B = L LT   мы получаем

  A z  =  l B z     = >     (L-1 A L-T) (LT z)  =  l  (LT z)  .

При этом собственные значения уравнения  A z = l B z  являются также собственными значениями уравнения  C y = l y,  где  симметричная матрица  C = L-1 A L-T ,  а   y = LT z .

В случае эрмитовых матриц  C = L-1 A L-H ,  а  y = LH z .

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

3.3. Восстановление собственных векторов для обобщенной проблемы собственных значений  z из полученных собственных векторов для линейной проблемы собственных значений  y может быть выполнено просто использованием соответствующих подпрограмм из пакета BLAS[20,21] уровня 2 или 3.

В таблице ниже показано, как каждая из трех видов проблем собственных значений может быть сведена к стандартной форме  C y = l y,  и как собственные векторы  z  исходной проблемы могут быть получены из собственных векторов  y  стандартной (линейной) проблемы.

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

  ----------------------------------------------------------------------------------------------------
                                                                                                      Восстановление
          Вид                                                                                            исходных
      проблемы            Разложение               Сведение                   собственных
                                                                                                             векторов
  -----------------------------------------------------------------------------------------------------
      Az = l Bz               B = LLT                  C = L-1 AL-T                    z = L-T y
                                    B = UT U                 C = U-T AU-1                   z = U-1 y
  -----------------------------------------------------------------------------------------------------
      ABz = l z               B = LLT                 C = LT AL                         z = L-T y
                                    B = UT U                C = UA UT                        z = U-1 y
  -----------------------------------------------------------------------------------------------------
      BAz = l z              B = LLT                  C = LT AL                         z = L y
                                   B = UT U                 C = UA UT                        z = UT y
  -----------------------------------------------------------------------------------------------------

Разложение Холецкого матрицы B выполняется с помощью базовой подпрограммы PDPOTRF.

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

Далее используется, необходимая для конкретного случая, одна из целевых подпрограмм Комплекса для решения линейной проблемы собственных значений ( PDSYEV1, PZHEGV1, PDSYEV2, PZHEGV2, PDSYEV3, PZHEGV3, PDSYEV4, PZHEGV4, PDSYEV5, PZHEGV5, PDSYEV6, PZHEGV6 ); см. п.2.

При необходимости вычисления собственных векторов, восстановление собственных векторов для обобщенной проблемы собственных значений по собственным векторам линейной проблемы собственных значений производится с помощью базовых подпрограмм PDTRSM ( для уравнений вида (1) или (2) ) или PDTRMM ( для уравнений вида (3) ).

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

PDSYGV1 -
PZHEGV1    		Вычисление всех собственных значений в обобщенной проблеме собственных значений вещественной симметричной (эрмитовой) матрицы с использованием  QL или QR - алгоритма
PDSYGV2 -
PZHEGV2    		Вычисление всех собственных значений и собственных векторов в обобщенной проблеме собственных значений вещественной симметричной (эрмитовой) матрицы с использованием  QL или QR - алгоритма
PDSYGV3 -
PZHEGV3    		Вычисление собственных значений в обобщенной проблеме собственных значений вещественной симметричной (эрмитовой) матрицы, принадлежащих заданному интервалу с использованием  QL или QR - алгоритма
PDSYGV4 -
PZHEGV4    		Вычисление собственных значений в обобщенной проблеме собственных значений, принадлежащих заданному интервалу, и соответствующих собственных векторов вещественной симметричной (эрмитовой) матрицы с использованием  QL или QR - алгоритма
PDSYGV5 -
PZHEGV5    		Вычисление собственных значений в обобщенной проблеме собственных значений, принадлежащих заданному интервалу индексов, вещественной симметричной (эрмитовой) матрицы с использованием  QL или QR - алгоритма
PDSYGV6 -
PZHEGV6    		Вычисление собственных значений в обобщенной проблеме собственных значений вещественной симметричной (эрмитовой) матрицы, принадлежащих заданному интервалу индексов, и соответствующих собственных векторов с использованием  QL или QR - алгоритма

4. Краткое описание содержания раздела Комплекса
для сингулярного разложения матриц

Пусть  A - вещественная матрица общего вида размера  m на n. Сингулярным разложением ( SVD ) матрицы  A  называется разложение вида 

                A  =  U S VT ,

где           U и V - ортогональные матрицы порядка  m и n  соответственно, а
                      S  =  diag (s1, ..., sr )

                 r  =  min (m, n) ,   s1 і ... і sr і 0.

Если  A - комплексная матрица, ее сингулярное разложение имеет вид

               A  =  U S VH ,

где         U и V - унитарные матрицы порядка  m и n  соответственно, а  S, также как и выше,
                           имеет вещественные диагональные элементы.

si называются сингулярными числами, первые  r  столбцов матрицы  V - правыми сингулярными векторами, а первые  r  столбцов матрицы  U - левыми сингулярными векторами матрицы  A.

Сингулярные числа и сингулярные векторы удовлетворяют условиям 

            A vi  =  si ui      и      AT ui  =  si vi   (или   AH ui  =  si vi ),

где  ui  и  vi   являются   i - ми столбцами матриц  U  и  V  соответственно.

Подробнее, вычисление состоит из следующих этапов.

4.1. Матрица  A  приводится к двухдиагональному виду A = U1 B V1T, если A  является вещественной (и A = U1 B V1H - если A - комплексная), где U1 и V1 являются ортогональными матрицами (или унитарными, если  A - комплексная), а  B  является вещественной и верхней двухдиагональной при  m і n (или нижней двухдиагональной при  m < n ).

4.2. Выполняется сингулярное разложение ( SVD ) двухдиагональной матрицы  B: 

         B  =  U2 S V2T ,

где     U2 и V2 - ортогональные матрицы, а S - диагональная (как указывалось выше).

Сингулярные числа матрицы  A совпадают с сингулярными числами матрицы  B. Сингулярными векторами матрицы  A  являются  U = U1 U2  и  V = V1 V2.

Приведение вещественной матрицы общего вида к двухдиагональной форме выполняется базовой подпрограммой  PDGEBRD, а сингулярное разложение матрицы  B  выполняется с использованием базовой подпрограммы  DBDSQR из пакета LAPACK [12].
Если требуется вычислить только сингулярные числа, то делается обращение к подпрограмме DLASQ1, которая выполняет это вычисление, используя  qd- алгоритм (см.[24]).

Базовая подпрограмма  PDGEBRD представляет  U1 и  V1 в виде произведений элементарных матриц отражения.
Для вещественной матрицы  A  матрицы  U1 и  V1 могут быть умножены на другие матрицы без явного формирования U1  и  V1  посредством использования базовой подпрограммы PDORMBR.

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

Если  m >> n, может быть более эффективно сначала выполнить QR-разложение матрицы  A, используя базовую подпрограмму  PDGEQRF, а затем выполнить сингулярное разложение квадратной матрицы  R порядка  n, т.к. если  A = Q R  и R = U S VT, то сингулярное разложение  A  задается формулой A = ( Q U ) S VT. Аналогично, если  m << n, может быть более эффективно сначала выполнить  LQ факторизацию  A, используя базовую подпрограмму  DGELQF. Эти предварительные факторизации выполняются самими целевыми подпрограммами для прямоугольных матриц PDGESVD4, PDGESVD5, PDGESVD6.

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

PDGESVD1 - Вычисление сингулярных чисел квадратной матрицы
PDGESVD2 - Сингулярное разложение квадратной матрицы с выдачей всех сингулярных чисел и сингулярных векторов
PDGESVD3 - Сингулярное разложение квадратной матрицы с выдачей сингулярных чисел и правых и/или левых сингулярных векторов
PDGESVD4 - Вычисление сингулярных чисел прямоугольной матрицы
PDGESVD5 - Сингулярное разложение прямоугольной матрицы с выдачей всех сингулярных чисел и сингулярных векторов
PDGESVD6 - Сингулярное разложение прямоугольной матрицы с выдачей сингулярных чисел и правых и/или левых сингулярных векторов

В Приложении 3_4 приводится конкретный пример обращения к целевой подпрограмме вычисления сингулярных чисел квадратной матрицы.

5. Разбор примера по использованию подпрограммы
для решения линейной симметричной проблемы собственных значений

Рассмотрим подробнее пример по использованию подпрограммы PDSYEV2, которая предназначена для вычисления всех собственных значений и собственных векторов симметричной матрицы общего вида  A. Данный пример приводится в описании этой подпрограммы в разделе "Пример использования" (в Приложении 3_1).

Пример демонстрирует решение проблемы собственных значений для симметричной матрицы  A размера 7 * 7. Ниже приводится ее изображение в символическом и числовом виде.

a11   a12   a13   a14   a15   a16   a17  1    0    0    0    0    0    1
a21   a22   a23   a24   a25   a26   a27  0    1    0    0    0    0    2
a31   a32   a33   a34   a35   a36   a37  0    0    1    0    0    0    3
a41   a42   a43   a44   a45   a46   a47  = 0    0    0    1    0    0    4
a51   a52   a53   a54   a55   a56   a57  0    0    0    0    1    0    5
a61   a62   a63   a64   a65   a66   a67  0    0    0    0    0    1    6
a71   a72   a73   a74   a75   a76   a77  1    2    3    4    5    6    7

Порядок матрицы N = 7.

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

Как говорилось в пп.4, 6.1 в первой части Методического пособия, базовые подпрограммы Комплекса реализуют специальные блочные алгоритмы обработки матриц, позволяющие повысить скорость решения задачи. При этом распределение плотных матриц по параллельных процессам происходит блочно - циклическим образом (см. пп.6.1, 6.2 из первой части пособия и п.4 во второй части пособия ).

Пусть матрица разбивается на блоки размером  MB * NB, где  MB = NB = 2.

1    0
0    1 		0    0
0    0 		0    0
0    0 		1
2
0    0
0    0 		1    0
0    1 		0    0
0    0 		3
4
0    0
0    0 		0    0
0    0 		1    0
0    1 		5
6
1    2 3    4 5    6 7

Пусть пропуск программы осуществляется на 4 процессах, которые образуют решетку 2 на 2, т.е. NPROW = 2 и NPCOL = 2 (см. пп.3, 4 в первой части пособия ).

После распределения полученных блоков по процессам решетки (по принятым правилам, которые описаны в первой части пособия в пп.6.1, 6.2 ), имеем следующую картинку распределения элементов матрицы  A по процессам (в символическом и числовом виде).
(Вверху указаны номера столбцов решетки процессов, слева - номера строк решетки процессов.)

0 1
0 a11  a12
a21  a22 		a15  a16
a25  a26 		a13  a14
a23  a24 		a17
a27
a51  a52
a61  a62 		a55  a56
a65  a66 		a53  a54
a63  a64 		a57
a67
1 a31  a32
a41  a42 		a35  a36
a45  a46 		a33  a34
a43  a44 		a37
a47
a71  a72 a75  a76 a73  a74 a77



0 1
0 1    0
0    1 		0    0
0    0 		0    0
0    0 		1
2
0    0
0    0 		1    0
0    1 		0    0
0    0 		5
6
1 0    0
0    0 		0    0
0    0 		1    0
0    1 		3
4
1    2 5    6 3    4 7

Из последнего рисунка видно, что процесс с координатами (0, 0) содержит локальный массив (часть матрицы  A) размера (4, 4), процесс (0, 1) - локальный массив размера (4, 3), процесс (1, 0) - локальный массив размера (3, 4), процесс (1, 1) - локальный массив размера (3, 3).

В Приложении 3_1 (раздел "Пример использования") приведен способ инициализации элементов распределенной матрицы  A с помощью обычных операторов присваивания (см. подпрограмму MATINIT), каждый из которых выполняется только на одном из параллельных процессов.

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

Такой способ инициализации распределенной матрицы  A явно не годится для матриц больших размеров. Поэтому в Приложении 3_2 приводится более простой способ инициализации матрицы  A с помощью обращения к служебной подпрограмме PDELSET, которая сама определяет (в соответствии с принятыми правилами блочно - циклического распределения), на какой из параллельных процессов и в какое место локальной памяти (для локальной части матрицы  A) следует поставить конкретный элемент глобальной матрицы  A с глобальными индексами   I, J  т.е.  A ( I, J ).

Оба приведенных в двух Приложениях способа распределения матрицы  A приводят к одному и тому же результату.

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

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

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

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

Полный текст, описанного в Приложении 3_3 примера использования подпрограммы PDSYEV2, можно получить также в системе Интернет по адресу: http://srcc.msu.su/num_anal/par_prog/ в разделе "Систематический каталог" (и здесь tdsyev25.zip).

В указанном zip - файле содержатся как файлы с головной программой, исходной матрицей и результатами, так и файлы со служебными подпрограммами PDLAREAD, PDLAWRITE, NUMROC.

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

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

Рассматривается решение проблемы собственных значений для эрмитовой матрицы  A размера 7 * 7.

Для простоты, демонстрация делается на примере матрицы, аналогичной рассмотренной в предыдущем пункте (п.5), т.е. подпрограмме PZHEEV4 подается частный случай эрмитовой матрицы, все мнимые части элементов которой равны 0.

Ниже приводится ее изображение в числовом виде.

1 0 0 0 0 0 1+i*0
0 1 0 0 0 0 2+i*0
0 0 1 0 0 0 3+i*0
0 0 0 1 0 0 4+i*0
0 0 0 0 1 0 5+i*0
0 0 0 0 0 1 6+i*0
1-i*0 2-i*0 3-i*0 4-i*0 5-i*0 6-i*0 7

Порядок матрицы N = 7.
Пусть интервал, в котором ищутся собственные значения ( VL, VU ) равен ( 2, 20 ).

Сначала зададим решетку процессов и распределим по ней исходную глобальную матрицу  A.

Пусть матрица разбивается на блоки размером  MB * NB, где  MB = NB = 2.

1      0
0      1 		0      0
0      0 		0      0
0      0 		1 + i*0
2 + i*0
0      0
0      0 		1      0
0      1 		0      0
0      0 		3 + i*0
4 + i*0
0      0
0      0 		0      0
0      0 		1      0
0      1 		5 + i*0
6 + i*0
1-i*0   2-i*0 3-i*0   4-i*0 5-i*0   6-i*0 7

Пусть пропуск программы осуществляется на 6 процессах, которые образуют решетку 3 на 2, т.е. NPROW = 3 и NPCOL = 2 (см. пп.3, 4 в первой части пособия ).

После распределения полученных блоков по процессам решетки (по принятым правилам, которые описаны в п.6.1 в первой части пособия ), имеем следующую картинку распределения элементов матрицы  A по процессам (в числовом виде).

(Вверху указаны номера столбцов решетки процессов, слева - номера строк решетки процессов.)

0 1 2
0 1     0
0     1 		 1+i*0
 2+i*0 		0     0
0     0 		0     0
0     0
0     0
0     0 		 5+i*0
 6+i*0 		0     0
0     0 		1     0
0     1
1 0     0
0     0 		 3+i*0
 4+i*0 		1     0
0     1 		0     0
0     0
1-i*0  2-i*0  7 3-i*0  4-i*0 5-i*0  6-i*0

Из последнего рисунка видно, что процесс с координатами (0, 0) содержит локальный массив (часть матрицы  A) размера (3, 4), процесс (0, 1) - локальный массив размера (2, 4), процесс (0, 2) - локальный массив размера (2, 4), процесс (1, 0) - локальный массив размера (3, 3), процесс (1, 1) - локальный массив размера (2, 3), процесс (1, 2) - локальный массив размера (2, 3).

Приведем здесь фортранные операторы второго, из рассмотренных в п.6.3 (в первой части пособия ), более простого способа реализации распределения на ФОРТРАНе исходной матрицы в локальную память каждого из параллельных процессов с помощью служебной подпрограммы PZELSET (предназначенной для комплексных матриц). 

*
      DO 2 J = 1, N
      DO 1 I = 1, N
      CALL PZELSET( A, I, J, DESCA, (0.0D0, 0.0D0) )
   1 CONTINUE
   2 CONTINUE
      DO 4 J = 1, N
      DO 3 I = 1, N
         IF( I.EQ.J ) THEN
            CALL PZELSET( A, I, J, DESCA, 1.0D0 )
         END IF
   3 CONTINUE
   4 CONTINUE
      J = N
      DO 5 I = 1, N
      CALL PZELSET( A, I, J, DESCA, DCMPLX( DBLE(I), 0.0D0) )
   5 CONTINUE
      I = N         
      DO 6 J = 1, N
      CALL PZELSET( A, I, J, DESCA, DCMPLX( DBLE(J), 0.0D0) )
   6 CONTINUE

Полностью весь фортранный текст рассматриваемого примера приведен в tzheev46.zip.

В результате получено одно собственное значение и соответствующий собственный вектор. 

  собственное число:

  W(1) = 14.0D0

  собственный вектор:

 Z( 1, 1 ) =  0.620173672946042268D-01 + i * 0.0D+0
 Z( 2, 1 ) =  0.124034734589208454D+00 + i * 0.0D+0
 Z( 3, 1 ) =  0.186052101883812687D+00 + i * 0.0D+0
 Z( 4, 1 ) =  0.248069469178416907D+00 + i * 0.0D+0
 Z( 5, 1 ) =  0.310086836473021155D+00 + i * 0.0D+0
 Z( 6, 1 ) =  0.372104203767625430D+00 + i * 0.0D+0
 Z( 7, 1 ) =  0.806225774829855024D+00 + i * 0.0D+0

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

7.Обработка ошибок и выдача диагностических сообщений.

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

7.1. Первый уровень представляют ситуации, диагностируемые самими подпрограммами Комплекса, как целевыми, так и базовыми, в частности, подпрограммами (используемого Комплексом) пакета PBLAS (см.).

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

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

Более подробно, в случае обнаружения ошибки в фактическом параметре, параметру INFO присваивается:

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

При диагностировании ошибки во входном параметре вызывается служебная подпрограмма обработки ошибок PXERBLA( ), которая выдает сообщение, содержащее имя подпрограммы, обнаружевшей ошибку, и номер ошибочного параметра, например:

"On entry to PDGESV parameter number 4 had an illegal value".

Обычная версия подпрограммы PXERBLA ( ), сделав выдачу диагностического сообщения, не останавливает выполнение программы. Это делается  из - за того, что некоторые "ошибки" могут быть исправимыми, и пользователю предоставляется возможность продолжить выполнение программы. В случае, если это не устраивает пользователя, он может использовать свой вариант PXERBLA ( ), добавив туда вызов подпрограммы BLACS_ABORT ( ) из пакета BLACS[17,18], которая прекратит выполнение программы пользователя.

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

Если же ошибка обнаружена в подпрограмме низкого уровня, она считается неисправимой, PXERBLA ( ) выдает сообщение и выполнение прерывается посредством вызова BLACS_ABORT ( ).

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

Подпрограммы пакета PBLAS[15,16] в целях эффективности выполняют проверку правильности только локальных параметров из своего списка параметров. Если ошибка обнаруживается хотя бы на одном из параллельных процессов текущего контекста, выполнение программы прерывается.

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

Отсутствие глобальной проверки в подпрограммах PBLAS'а является следствием высокой стоимости (в смысле ухудшения производительности) таких проверок.

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

Это бывает по следующим основным причинам:

Например, если используется подпрограмма (PSGESV2) для решения системы уравнений с вещественной матрицей одинарной точности, которая является близкой к вырожденной, может быть зафиксирована явная вырожденность на  i - ой  итерации  LU - факторизации.

В этом случае либо INFO присваивается значение  i,  либо (что более вероятно) может быть вычислена оценка обратного числа обусловленности, которая окажется меньше, чем относительная машинная точность; в этом случае INFO полагается равным  (n + 1).

Если возникла ситуация, когда выдается INFO > 0, управление всегда возвращается в вызывающую подпрограмму, а подпрограмма PXERBLA ( ) не вызывается, и сообщение об ошибке не выдается.

Поэтому при выходе из подпрограммы Комплекса всегда следует проверять значение параметра INFO.

Ошибка с INFO > 0 может получиться, например, в одном из следующих случаев.

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

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

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

Этот уровень устанавливается посредством использования макроса препроцессора языка Си  BlacsDebugLvl.

Чтобы узнать, какой уровень отладки используется в BLACS'е в данном случае, можно использовать подпрограмму  BLACS_GET ( ).

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

Но большинство параметров BLACS не проверяет,  из - за того, что это приведет к ухудшению производительности.

Поэтому пользователям рекомендуется при связывании своей программы с библиотекой BLACS'а устанавливать при компиляции отладочный уровень  1, до тех пор пока идет процесс отладки программы. При этом проверяется большая часть параметров. Кроме того, выдаются и некоторые другие полезные сообщения. Например, пользователь будет предупрежден, если какой - либо процесс посылает сообщение сам себе. Такая ситуация допустима, однако, свидетельствует о плохом программировании и требует достаточно большого буферного пространства, потому что в этом случае может произойти "зависание" при передаче больших сообщений.

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

Все сообщения BLACS'а можно разделить на три следующих вида.

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

Сообщение выдается в виде: 

BLACS WARNING  'текст сообщения'
from {‹p›,‹q›},  pnum = ‹pnum›,  Contxt = ‹ictxt›,  on line ‹#›  of file  '‹fname›', 
здесь
{‹p›, ‹q›} - координаты процесса, выдавшего сообщение, в решетке процессов;
‹pnum› - номер процесса, выдаваемый как выходное значение первого параметра подпрограммой BLACS_PINFO;
‹ictxt› - целое значение контекста. Обратите внимание на то, что это значение не является одинаковым для всех процессов решетки. Например, процесс {0, 0} может иметь ictxt = 0, а процесс {0, 1} имеет ictxt = 1 для того же самого контекста. Однако, ‹pnum› и ‹ictxt› вместе обеспечивают однозначную идентификацию процесс/контекст;
‹#› - номер строки в файле с именем ‹fname›, который выдает сообщение;
‹fname› - имя файла, в котором содержится подпрограмма, выдавшая сообщение.

б). Ошибки. BLACS обнаружил ошибку. Выдается сообщение об ошибке посредством внутренней подпрограммы  BlacsErr, и выполнение программы останавливается посредством вызова подпрограммы  BLACS_ABORT ( ). Сообщение при этом выдается в таком же виде, как и в случае предупреждения, за исключением того, что слова  BLACS WARNING заменяются на слова  BLACS ERROR.

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

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

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

в). Системная ошибка. BLACS получил сообщение об ошибке от операционной системы. Он передает его пользователю и завершает выполнение программы.

7.3. К третьему уровню ошибок, относятся системные ошибки.

При работе с подпрогрммами Комплекса иногда могут возникать сообщения об ошибках, передаваемые операционной системой. В предыдущем п.7.2_в) описана реакция BLACS'а на эту ситуацию.

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

Например, если используется PVM BLACS, будет выдан номер ошибки в стиле PVM. Тогда для понимания сообщения об ошибке, т.е. перевода номера сообщения в соответствующий текст, может потребоваться, например, краткое руководство по PVM.

ЛИТЕРАТУРА

  1. Антонов А.С. Параллельное программирование с использованием технологии MPI: Учебное пособие: -М.: Изд-во МГУ, 2004.-71с.
  2. Антонов А.С. Практический курс MPI. MPI: A Message-Passing Interface Standard (Version 1.1).
  3. Библиотека численного анализа НИВЦ МГУ, http://www.srcc.msu.su/num_anal/
  4. Воеводин В.В. Математические модели и методы в параллельных процессах. -М.: Наука, 1986.-296с.
  5. Воеводин В.В. Математические основы параллельных вычислений. -М.: Изд-во МГУ, 1991.-345с.
  6. Воеводин В.В. Параллельные структуры алгоритмов и программ. -М.: ОВМ АН СССР, 1987.-148с.
  7. Воеводин В.В., Воеводин Вл.В. Параллельные вычисления. -С.-П. -"БХВ-Петербург", 2002.-608с.
  8. Дацюк В.Н., Букатов А.А., Жегуло А.И. Введение в организацию и методы программирования многопроцессорных вычислительных систем (методическое пособие, часть I) -Ростов-на-Дону: Изд-во РГУ, 2000.-36с.
  9. Дацюк В.Н., Букатов А.А., Жегуло А.И. Среда параллельного программирования MPI (методическое пособие, часть II) -Ростов-на-Дону: Изд-во РГУ, 2000.-65с.
  10. Немнюгин С.А., Стесик О.Л. "Параллельное программирование для многопроцессорных вычислительных систем", "БХВ-Петербург", 2002.-400с.
  11. Корнеев В.Д. "Параллельное программирование в MPI", Новосибирск: Изд-во СО РАН, 2000.-213с.
  12. E. Anderson, Z. Bai, C. Bischof, J.Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. Mckenney, S. Ostrouchov, and D. Sorensen, LAPACK Users' Guide, Society for Industrial and Applied Mathematics, Philadelphia, PA, second ed., 1995.
  13. L. S. Blackford, J. Choi, A. Cleary, J. Demmej, I. Dhillon, J. J. Dongarra, S. Hammarling, G. Henry, A. Petitet, K. Stanley, D. W. Walker, and R. C. Whaley, ScaLAPACK: A portable linear algebra library for distributed memory computers - design issues and performance, in Proceedings of Supercomputing '96, Sponsored by ACM SIGARCH and IEEE Computer Society, 1996. (ACM Order Number: 415962, IEEE Computer Society Press Order Number: RS00126. http://www.supercomp.org/sc96/proceedings/).
  14. http://www.netlib.org/scalapack/index.html/.
  15. J. Choi, J. Dongarra, and D. Walker, PB-BLAS: A Set of Parallel Block Basic Linear Algebra Subroutines, Practice and Experience, 8 (1996), pp. 517-535.
  16. http://www.netlib.org/scalapack/html/pblas_qref.html/.
  17. J. Dongarra, R. van de Geun, and R.C. Whaley, Two dimensional basic linear algebra communication subprograms, in Environments and Tools for Parallel Scientific Computing, Advances in Parallel Computing, J. Dongarra and B. Tourancheau, eds., vol. 6, Elsevier Science Publishers B.V., 1993, pp. 31-40.
  18. http://www.netlib.org/blacs/
  19. G. Golub and C.F. Van Loan, Matrix Computations, Johns Hopkins University Press, Baltimore, MD, third ed., 1996.
  20. C.L. Lawson, R.J. Hanson, D. Kincaid, and F.T. Krogh, Basic linear algebra subprograms for Fortran usage, ACM Trans. Math. Soft., 5 (1979), pp. 308-323.
  21. http://www.netlib.org/blas/
  22. W. Lichtenstein and S.L. Johnsson, Block-cyclic dense linear algebra, SIAM J. Sci. Stat. Comput. 14 (1993), pp. 1259-1288.
  23. L. Prylli and B. Tourancheau, Efficient block cyclic data redistribution, in EUROPAR'96, vol. 1 of Lecture Notes in Computer Science, Springer-Verlag, 1996, pp. 155-165.
  24. Fernando K.V. and Parlett B.N. Accurate singular values and differential qd algorithms// Numer.Math. 1994, Vol-67, No. 2, pp. 191-230.
  25. Форсайт Дж., Малькольм М., Моулер К. Машинные методы математических вычислений. М.: Изд-во МИР, 1980.

Приложение 3_1

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

Подпрограмма:  PDSYEV2

Назначение

Вычисление всех собственных значений и собственных векторов вещественной симметричной матрицы

Математическое описание

На первом этапе подпрограмма приводит симметричную матрицу  A к трехдиагональной форме методом отражений:
A =  Q*D*QT, где Q - ортогональная матрица, а D - симметричная трехдиагональная матрица. На втором этапе неявным  QL или  QR алгоритмом находятся собственные значения матрицы  D, которые совпадают с собственными значениями матрицы A, поскольку матрицы  A и  D ортогонально подобны. Собственные значения и собственные векторы вычисляются с машинной точностью.

Литература:
http://www.srcc.msu.su/num_anal/par_prog/

Использование 

         CALL  PDSYEV2 ( UPLO, N, A, IA, JA, DESCA, W,
                                         Z, IZ, JZ, DESCZ, WORK, LWORK, INFO)

   Параметры
UPLO - если UPLO = ' U ' - задается верхний треугольник матрицы  A;
если UPLO = ' L ' - задается нижний треугольник матрицы  A
(глобальный входной параметр, тип символьный);
N - порядок распределенной подматрицы sub (A) (см. п.5.3); N і 0 (глобальный входной параметр, тип целый);
A - массив двойной точности, распределенный по процессам блочно - циклическим образом (см. пп.5.1, 6.1, 6.2), глобальная размерность которого (N, N), а локальная размерность ( LLD_A, LOCc ( JA + N - 1));
на входе - это симметричная матрица  A;
если UPLO = ' U ', то используется только верхняя треугольная часть матрицы  A;
если UPLO = ' L ', то используется только нижняя треугольная часть матрицы  A;
на выходе - нижний треугольник (при UPLO = ' L ' ) или верхний треугольник (при UPLO = ' U ' ) матрицы  A, включая диагональ, не сохраняется (локальный входной параметр, локальное рабочее пространство);
IA - глобальный номер строки матрицы  A, который указывает на начало подматрицы (см. п.5.3) (глобальный входной параметр, тип целый);
JA - глобальный номер столбца матрицы  A, который указывает на начало подматрицы (см. п.5.3) (глобальный входной параметр, тип целый);
DESCA - дескриптор распределенной матрицы  A (одномерный массив длины DLEN);
если тип дескриптора одномерный ( DTYPE_A = 501),  DLEN і 7,
если тип дескриптора двумерный    ( DTYPE_A = 1),      DLEN і 9;
DESCA содержит информацию о размещении  A в памяти; полное описание DESCA см. в п.5.2 (глобальный и локальный входной параметр, тип целый);
W - одномерный массив двойной точности длины N; если  INFO = 0, собственные значения располагаются в нем по возрастанию (глобальный выходной параметр);
Z - массив двойной точности с глобальной размерностью ( N, N ) и с локальной размерностью ( LLD_Z, LOCc ( JZ + N - 1)); содержит ортонормированные собственные векторы симметричной матрицы  A (локальный выходной параметр);
IZ - глобальный номер строки матрицы Z, который указывает на начало подматрицы sub (Z) (глобальный входной параметр, тип целый);
JZ - глобальный номер столбца матрицы Z, который указывает на начало подматрицы sub (Z) (глобальный входной параметр, тип целый);
DESCZ - дескриптор распределенной матрицы  Z (одномерный массив длины DLEN);
если тип дескриптора одномерный ( DTYPE_Z = 501),  DLEN і 7,
если тип дескриптора двумерный    ( DTYPE_Z = 1),      DLEN і 9;
DESCZ содержит информацию о размещении  Z в памяти; полное описание DESCZ см. в п.5.2; DESCZ (CTXT_) должен быть равен DESCA (CTXT_)
(глобальный и локальный входной параметр, тип целый);
WORK - одномерный рабочий массив двойной точности длины LWORK;
на выходе, в элементе WORK (1) возвращается необходимая длина рабочего пространства
(локальное рабочее пространство, локальный выходной параметр);
LWORK - задаваемая длина рабочего пространства WORK;
Минимальная требуемая величина LWORK вычисляется самой подпрограммой, если к ней обратиться со значением LWORK = -1; при этом подпрограмма не производит никаких других вычислений, а требуемое значение LWORK возвращается в элементе массива WORK (1);
(локальный или глобальный входной параметр, тип целый);
INFO - целая переменная, диагностирующая результат работы подпрограммы (глобальный выходной параметр)
= 0 - успешное завершение работы;
< 0 - если  i - ый фактический параметр подпрограммы является массивом и его  j - ый элемент имеет недопустимое значение, тогда INFO = - ( i * 100 + j ), если  i - ый фактический параметр является скаляром и имеет недопустимое значение, тогда INFO = - i;
> 0 - если INFO = i, где 1 Ј i Ј N, то  i - ое собственное значение не может быть вычислено с машинной точностью после выполнения 30 * N итераций;
если INFO = N+1, то подпрограмма не гарантирует точности вычисленных собственных значений. 
   Версии
PSSYEV2 -   PCHEEV2     PZHEEV2     вычисление всех собственных значений и собственных векторов симметричной (эрмитовой) матрицы для случаев вещественных данных одинарной точности, комплексных данных одинарной точности, комплексных данных двойной точности соответственно 
   Вызываемые подпрограммы

Здесь указаны только базовые подпрограммы (2 - ого уровня), которые вызываются из целевой подпрограммы (1 - ого уровня).

PDSYTRD -
PZHETRD    		приведение симметричной (эрмитовой) матрицы к трехдиагональной форме методом отражений
DSTEQR2 -
ZSTEQR2    		вычисление всех собственных значений и, возможно, собственных векторов симметричной трехдиагональной матрицы
PDORMTR -
PZUNMTR    		умножение матрицы общего вида на матрицу отражения, вычисленную подпрограммой PDSYTRD (PZHETRD)
PDLAMCH - вычисление машинных параметров для арифметики с плавающей запятой
PDLASET - внедиагональные элементы матрицы полагаются равными ALFA, а диагональные элементы - равными BETA
PDLANSY - вычисляет значения 1 - ой нормы, нормы Фробениуса, бесконечной нормы или наибольшее абсолютное значение любого элемента вещественной симметричной, комплексной симметричной или эрмитовой матрицы
PDLASCL - умножает вещественную матрицу на вещественный скаляр 
   Замечания по использованию
  1.  В подпрограммах PSSYEV2, PCHEEV2, PZHEEV2 параметры A, WORK и Z имеют тип REAL, COMPLEX и DOUBLE COMPLEX соответственно, а параметр W - тип REAL, REAL и DOUBLE PRECISION соответственно
Список параметров подпрограммы PZHEEV2 имеет следующий вид:
PZHEEV2(UPLO, N, A, IA, JA, DESCA, W, Z, IZ, JZ, DESCZ, WORK, LWORK, RWORK, LRWORK, INFO),
где дополнительные параметры RWORK и LRWORK означают:
RWORK - одномерный рабочий массив двойной точности длины LRWORK; на выходе в элементе RWORK(1) возвращается необходимая длина рабочего массива RWORK (локальное рабочее пространство, локальный выходной параметр);
LRWORK - задаваемая длина рабочего пространства RWORK; LRWORK і 2 * N + 2 * N - 2 (локальный входной параметр, тип целый);
  2.  Используются подпрограммы BLACS_GRIDINFO ( из пакета BLACS),
DCOPY, DSCAL ( из пакета BLAS),
LSAME, INDXG2P, NUMROC, CHK1MAT, PCHK2MAT, PXERBLA, SL_GRIDRESHAPE, PDELGET ( из библиотеки ScaLAPACK_TOOLS)

Пример использования

Вычисление всех собственных значений и собственных векторов симметричной матрицы  A, которая имеет вид:

 1  0  0  0  0  0  1
 0  1  0  0  0  0  2
 0  0  1  0  0  0  3
 0  0  0  1  0  0  4
 0  0  0  0  1  0  5
 0  0  0  0  0  1  6
 1  2  3  4  5  6  7

Порядок матрицы N = 7.
Пусть матрица разбивается на блоки с размером блоков NB = 2.

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

Фрагмент фортранного текста вызывающей программы.
( принятые в тестах обозначения можно посмотреть в первой части методического пособия [п.8.3]).

В Приложении 3_2 приведен другой (более простой) способ инициализации и распределения матрицы  A посредством обращения к служебной подпрограмме п.PDELSET. 

      PROGRAM TDSYEV21
      include 'mpif.h'
      INTEGER                        LWORK, MAXN, LDA, MAXPROCS, NOUT
      DOUBLE PRECISION   ZERO, MONE
      PARAMETER                 ( MAXN = 100, LDA = 100, LWORK = 264,
     $                                      MAXPROCS = 512, NOUT =6, ZERO = 0.0D+0,
     $                                      MONE = -1.0D+0 )
      CHARACTER                UPLO
      PARAMETER                 ( UPLO = 'U' )

      INTEGER                        CTXT, I, IAM, INFO, MYCOL, MYROW, N, NB,
     $                                       NPCOL, NPROCS, NPROW, IA, JA, IZ, JZ

      INTEGER                        DESCA( 9 ), DESCZ( 9 )
      DOUBLE PRECISION   A( LDA, LDA ), Z( LDA, LDA ), W( MAXN ),
     $                                       WORK( LWORK ), WORK1( 7 )

      EXTERNAL                    BLACS_EXIT, BLACS_GET, BLACS_GRIDEXIT,
     $                                       BLACS_GRIDINFO, BLACS_GRIDINIT, BLACS_PINFO,
     $                                       BLACS_SETUP, DESCINIT, MATINIT, PDLAPRNT,
     $                                       PDSYEV2

      N = 7
      NB = 2
      NPROW = 2
      NPCOL = 2
      IA = 1
      JA = 1
      IZ = 1
      JZ = 1

      CALL BLACS_PINFO( IAM, NPROCS )
      IF( ( NPROCS.LT.1 ) ) THEN
         CALL BLACS_SETUP( IAM, NPROW*NPCOL )
      END IF

      CALL BLACS_GET( -1, 0, CTXT )
      CALL BLACS_GRIDINIT( CTXT, 'R', NPROW, NPCOL )
      CALL BLACS_GRIDINFO( CTXT, NPROW, NPCOL, MYROW, MYCOL )

      IF( MYROW.EQ.-1 ) GO TO 20

      CALL DESCINIT( DESCA, N, N, NB, NB, 0, 0, CTXT, LDA, INFO )
      CALL DESCINIT( DESCZ, N, N, NB, NB, 0, 0, CTXT, LDA, INFO )
*
* Построение матрицы A
*
      CALL MATINIT( N, A, IA, JA, DESCA )
*
*  Вычисление всех собственных значений и собственных векторов
*
      CALL PDSYEV2( UPLO, N, A, IA, JA, DESCA, W,
     $                              Z, IZ, JZ, DESCZ, WORK, LWORK, INFO)
*
      CALL BLACS_GRIDEXIT( CTXT )
 20 CONTINUE

      CALL BLACS_EXIT( 0 )
      STOP
      END
*
      SUBROUTINE MATINIT( N, A, IA, JA, DESCA )
*
* MATINIT генерирует и распределяет матрицу A по решетке процессов
*
      INTEGER                         LLD_, CTXT_
      PARAMETER                  ( LLD_ = 9, CTXT_ = 2 )
      INTEGER                         IA, JA, N
      INTEGER                         DESCA( * )
      DOUBLE PRECISION    A( * )
      EXTERNAL                     BLACS_GRIDINFO
      INTEGER                         MXLLDA, MYCOL, MYROW, NPCOL, NPROW
      DOUBLE PRECISION    C, K
*
* Выполняемые операторы
*
      CALL BLACS_GRIDINFO( DESCA( CTXT_ ), NPROW, NPCOL, MYROW, MYCOL )
*
      C = 0.0D0
      K = 1.0D0
*
      MXLLDA = DESCA( LLD_ )
*
* Локальная часть матрицы A для процесса (0, 0)
*
      IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
         A( 1 ) = K
         A( 2 ) = C
         A( 3 ) = C
         A( 4 ) = C
         A( 1+MXLLDA ) = C
         A( 2+MXLLDA ) = K
         A( 3+MXLLDA ) = C
         A( 4+MXLLDA ) = C
*
         A( 1+2*MXLLDA ) = C
         A( 2+2*MXLLDA ) = C
         A( 3+2*MXLLDA ) = K
         A( 4+2*MXLLDA ) = C
*
         A( 1+3*MXLLDA ) = C
         A( 2+3*MXLLDA ) = C
         A( 3+3*MXLLDA ) = C
         A( 4+3*MXLLDA ) = K
*
* Локальная часть матрицы A для процесса (0, 1)
*
      ELSE IF( MYROW.EQ.0 .AND. MYCOL.EQ.1 ) THEN
         A( 1 ) = C
         A( 2 ) = C
         A( 3 ) = C
         A( 4 ) = C
*
         A( 1+MXLLDA ) = C
         A( 2+MXLLDA ) = C
         A( 3+MXLLDA ) = C
         A( 4+MXLLDA ) = C
*
         A( 1+2*MXLLDA ) = K
         A( 2+2*MXLLDA ) = 2.0D0
         A( 3+2*MXLLDA ) = 5.0D0
         A( 4+2*MXLLDA ) = 6.0D0
*
* Локальная часть матрицы A для процесса (1, 0)
*
      ELSE IF( MYROW.EQ.1 .AND. MYCOL.EQ.0 ) THEN
         A( 1 ) = C
         A( 2 ) = C
         A( 3 ) = K
*
         A( 1+MXLLDA ) = C
         A( 2+MXLLDA ) = C
         A( 3+MXLLDA ) = 2.0D0
*
         A( 1+2*MXLLDA ) = C
         A( 2+2*MXLLDA ) = C
         A( 3+2*MXLLDA ) = 5.0D0
*
         A( 1+3*MXLLDA ) = C
         A( 2+3*MXLLDA ) = C
         A( 3+3*MXLLDA ) = 6.0D0
*
* Локальная часть матрицы A для процесса (1, 1)
*
      ELSE IF( MYROW.EQ.1 .AND. MYCOL.EQ.1 ) THEN
         A( 1 ) = K
         A( 2 ) = C
         A( 3 ) = 3.0D0
*
         A( 1+MXLLDA ) = C
         A( 2+MXLLDA ) = K
         A( 3+MXLLDA ) = 4.0D0
*
         A( 1+2*MXLLDA ) = 3.0D0
         A( 2+2*MXLLDA ) = 4.0D0
         A( 3+2*MXLLDA ) = 7.0D0
*
      END IF
*
      RETURN
      END

Результаты:

 Значение  INFO  =   0

 Собственные значения:

  W(1)= -6.D+00
  W(2)=  1.D+00
  W(3)=  1.D+00
  W(4)=  1.D+00
  W(5)=  1.D+00
  W(6)=  1.D+00
  W(7)= 14.D+00

 Собственные векторы:

      Z(1,1)=  -0.845154254728516519D-01
      Z(2,1)=  -0.169030850945703026D+00
      Z(3,1)=  -0.253546276418554983D+00
      Z(4,1)=  -0.338061701891406607D+00
      Z(5,1)=  -0.422577127364258176D+00
      Z(6,1)=  -0.507092552837109856D+00
      Z(7,1)=   0.591607978309961591D+00

      Z(1,2)=  -0.255220953937281947D+00
      Z(2,2)=  -0.609603954124120162D+00
      Z(3,2)=   0.741527088114496191D+00
      Z(4,2)=   0.804104054667498236D-02
      Z(5,2)=  -0.108155528220763852D+00
      Z(6,2)=  -0.402564872068081891D-01
      Z(7,2)=   0.111022302462515654D-15

      Z(1,3)=   0.956284634899674191D+00
      Z(2,3)=  -0.109002334095028616D+00
      Z(3,3)=   0.218161506024083796D+00
      Z(4,3)=  -0.544996702163083949D-01
      Z(5,3)=  -0.112577050653521749D+00
      Z(6,3)=  -0.101980091774837717D+00
      Z(7,3)=   0.277555756156289135D-16

      Z(1,4)=   0.000000000000000000D+00
      Z(2,4)=  -0.986075661850500072D-02
      Z(3,4)=  -0.722934823299492685D-01
      Z(4,4)=   0.898666233992495123D+00
      Z(5,4)=  -0.239109239331969259D+00
      Z(6,4)=  -0.360419463180546074D+00
      Z(7,4)=   0.000000000000000000D+00

      Z(1,5)=   0.000000000000000000D+00
      Z(2,5)=  -0.729977775740844881D-02
      Z(3,5)=  -0.823564492566095790D-01
      Z(4,5)=   0.493401310222317258D-01
      Z(5,5)=  -0.759354818204546977D+00
      Z(6,5)=   0.643513745036408902D+00
      Z(7,5)=   0.000000000000000000D+00

      Z(1,6)=   0.969172365326708696D-01
      Z(2,6)=  -0.756567872068086622D+00
      Z(3,6)=  -0.540030650666323941D+00
      Z(4,6)=   0.105382651046173337D+00
      Z(5,6)=   0.259057033086912525D+00
      Z(6,6)=   0.219915781663869714D+00
      Z(7,6)=   0.222044604925031308D-15

      Z(1,7)=  -0.620173672946042337D-01
      Z(2,7)=  -0.124034734589208551D+00
      Z(3,7)=  -0.186052101883812715D+00
      Z(4,7)=  -0.248069469178416879D+00
      Z(5,7)=  -0.310086836473021099D+00
      Z(6,7)=  -0.372104203767625319D+00
      Z(7,7)=  -0.806225774829854913D+00

Приложение 3_2

Пример составления головной программы
для использования целевой подпрограммы из Приложения 3_1.

      PROGRAM TDSYEV21
*
*     Тест к подпрограмме   PDSYEV2
*
      include 'mpif.h'
* Именованные константы - параметры задачи
      INTEGER                        LWORK, MAXN, LDA, MAXPROCS, NOUT
      DOUBLE PRECISION   ZERO, MONE
      PARAMETER                 ( LDA = 100, LWORK = 70, MAXN = 100,
*      PARAMETER                 ( LDA = 100, LWORK = -1, MAXN = 100,
     $                                       MAXPROCS = 512, NOUT =6, ZERO = 0.0D+0,
     $                                       MONE = -1.0D+0 )
*
      CHARACTER          UPLO
      PARAMETER          ( UPLO = 'L' )
*
* Локальные переменные
      INTEGER            CTXT, I, IAM, INFO, MYCOL, MYROW, N, NB,
     $                            NPCOL, NPROCS, NPROW, IA, JA, IZ, JZ
*
* Локальные массивы
      INTEGER                       DESCA( 9 ), DESCZ( 9 )
      DOUBLE PRECISION   A( LDA, LDA ), W( MAXN ),
     $                                      WORK( LWORK ), Z( LDA, LDA ), WORK1( 7 )
*     $                                      WORK( 70 ), Z( LDA, LDA ), WORK1( 7 )
*
* Внешние подпрограммы
      EXTERNAL           BLACS_EXIT, BLACS_GET, BLACS_GRIDEXIT,
     $                              BLACS_GRIDINFO, BLACS_GRIDINIT, BLACS_PINFO,
     $                              BLACS_SETUP, DESCINIT, PDMATINIT, PDLAPRNT,
     $                              PDSYEV2
*
* Выполняемые операторы
*
* Постановка задачи
*
      N = 7
      NB = 2
      NPROW = 2
      NPCOL = 2
      IA = 1
      JA = 1
      IZ = 1
      JZ = 1
*
* Инициализация пакета BLACS
*
      CALL BLACS_PINFO( IAM, NPROCS )
      IF( ( NPROCS.LT.1 ) ) THEN
         CALL BLACS_SETUP( IAM, NPROW*NPCOL )
      END IF
*
* Инициализация контекста BLACS'а и решетки процессов
*
      CALL BLACS_GET( -1, 0, CTXT )
      CALL BLACS_GRIDINIT( CTXT, 'R', NPROW, NPCOL )
      CALL BLACS_GRIDINFO( CTXT, NPROW, NPCOL, MYROW, MYCOL )
*
* Если процесс не является частью данного контекста,
* переход на конец программы
*
      IF( MYROW.EQ.-1 ) GO TO 20
*
* Инициализация дескрипторов для матриц A и Z
*
      CALL DESCINIT( DESCA, N, N, NB, NB, 0, 0, CTXT, LDA, INFO )
      CALL DESCINIT( DESCZ, N, N, NB, NB, 0, 0, CTXT, LDA, INFO )
*
* Построение матрицы A
*
      CALL PDMATINIT( N, A, IA, JA, DESCA, INFO )
*
* Печать входных данных
*
      IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
         WRITE( NOUT, FMT = 92 )
         WRITE( NOUT, FMT = 98 )N, N, NB, NB
         WRITE( NOUT, FMT = 97 )NPROW*NPCOL, NPROW, NPCOL
         WRITE( NOUT, FMT = * ) ' Исходные данные '
         WRITE( NOUT, FMT = 88 ) UPLO
         WRITE( NOUT, FMT = * ) '  DESCA'
         WRITE( NOUT, FMT = 85 )DESCA
         WRITE( NOUT, FMT = * ) '  DESCZ'
         WRITE( NOUT, FMT = 85 )DESCZ
         WRITE( NOUT, FMT = * ) ' Матрица A'
      END IF
 85 FORMAT( / 9I5 /)
*  84 FORMAT( / ' WORK(1) = ', D16.5 /)
 88 FORMAT( / ' UPLO = ', A /)

      CALL PDLAPRNT( N, N, A, IA, JA, DESCA, 0, 0, 'A', NOUT, WORK1 )
*
* Вычисление собственных значений и собственных векторов
*
      CALL PDSYEV2( UPLO, N, A, IA, JA, DESCA, W, Z, IZ, JZ,
     $                             DESCZ, WORK, LWORK, INFO )
*
* Печать результатов
*
      IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
         WRITE( NOUT, FMT = * ) ' Результаты '
         WRITE( NOUT, FMT = 96 )INFO
*         WRITE( NOUT, FMT = 84 )WORK(1)
         WRITE( NOUT, FMT = * ) ' собственные числа'
         DO 10 I = 1, N
            PRINT *, ' W(', I, ')=', W( I ), ';'
 10    CONTINUE
         WRITE( NOUT, FMT = * ) ' собственные векторы'
      END IF
*
* Печать собственных векторов
*
      CALL PDLAPRNT( N, N, Z, IZ, JZ, DESCZ, 0, 0, 'Z', NOUT, WORK1 )
*
      CALL BLACS_GRIDEXIT( CTXT )
 20 CONTINUE

      CALL BLACS_EXIT( 0 )

 99 FORMAT( 'W=diag([', 4D16.12, //']);' )

 92 FORMAT( / 'Тест к подпрограмме PDSYEV2' )
 98 FORMAT( / ' Вычисление всех собственных значений и собственных',
     $  ' векторов матрицы  A,' //
     $      ' где A - симметричная матрица ',
     $      I2, ' на ', I2,', разбитая на блоки', I2, ' на ', I2, /)
 97 FORMAT( / ' Запуск на ', I2, ' процессах,',
     $       ' образующих решетку размером ', I2, ' на ', I2 /)
 96 FORMAT( / ' Значение INFO  = ', I6 /)
      STOP
      END
*
      SUBROUTINE PDMATINIT( N, A, IA, JA, DESCA, INFO )
*
* PDMATINIT генерирует и распределяет матрицу A по решетке процессов
*
      INTEGER                        IA, INFO, JA, N
      INTEGER                        DESCA( * )
      DOUBLE PRECISION   A( * )
      INTEGER                        I, J, MYCOL, MYROW, NPCOL, NPROW
      EXTERNAL                     PDELSET
      INTRINSIC                     DBLE
*
      INFO = 0
*
      IF( IA.NE.1 ) THEN
         INFO = -3
      ELSE IF( JA.NE.1 ) THEN
         INFO = -4
      END IF
*
      DO 2  J = 1, N
      DO 1  I = 1, N
      CALL PDELSET( A, I, J, DESCA, 0.0D0)
   1 CONTINUE
   2 CONTINUE
      DO 4  J = 1, N
      DO 3  I = 1, N
         IF( I.EQ.J ) THEN
            CALL PDELSET( A, I, J, DESCA, 1.0D0)
         END IF
   3 CONTINUE
   4 CONTINUE
      J = N
      DO 5  I = 1, N
      CALL PDELSET( A, I, J, DESCA, DBLE(I) )
   5 CONTINUE
      I = N
      DO 6  J = 1, N
      CALL PDELSET( A, I, J, DESCA, DBLE(J) )
   6 CONTINUE
*
      RETURN
      END

Приложение 3_3

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

Содержимое файла с исходной матрицей приведено следом за головной программой.

      PROGRAM TDSYEV25
*
*     Тест к подпрограмме   PDSYEV2
*
      include 'mpif.h'
* Именованные константы - параметры задачи ..
      INTEGER                    LWORK, NOUT, MEMSIZE, NIN, NPR
      PARAMETER             ( MEMSIZE = 20000, LWORK = 70,
     $                                   NOUT =13, NIN = 11, NPR = 6 )
*
      CHARACTER           UPLO
      PARAMETER          ( UPLO = 'L' )
*
* Локальные переменные
      INTEGER            CTXT, I, IAM, INFO, IPA, IPZ, IPW,
     $                            IPWORK, MYCOL, MYROW, NP, NQ,
     $                            N, NB, NPCOL, NPROCS, NPROW, IA, JA, IZ, JZ
*
* Локальные массивы
      INTEGER                         DESCA( 9 ), DESCZ( 9 )
      DOUBLE PRECISION    MEM( MEMSIZE ), WORK1(7)
*
* Внешние подпрограммы
      EXTERNAL           BLACS_EXIT, BLACS_GET, BLACS_GRIDEXIT,
     $                              BLACS_GRIDINFO, BLACS_GRIDINIT, BLACS_PINFO,
     $                              BLACS_SETUP, DESCINIT, PDLAPRNT, PDLAREAD,
     $                              PDLAWRITE, NUMROC, PDSYEV2
*
* Выполняемые операторы
*
* Постановка задачи
*
      N = 7
      NB = 2
      NPROW = 2
      NPCOL = 2
      IA = 1
      JA = 1
      IZ = 1
      JZ = 1
*
* Инициализация пакета BLACS
*
      CALL BLACS_PINFO( IAM, NPROCS )
*
        IF( NPROCS.LT.1 ) THEN
           CALL BLACS_SETUP( IAM, NPROW*NPCOL )
        END IF
*
* Инициализация контекста BLACS'а и решетки процессов
*
      CALL BLACS_GET( -1, 0, CTXT )
      CALL BLACS_GRIDINIT( CTXT, 'R', NPROW, NPCOL )
      CALL BLACS_GRIDINFO( CTXT, NPROW, NPCOL, MYROW, MYCOL )
*
* Если процесс не является частью данного контекста,
* переход на конец программы
*
      IF( MYROW.GE.NPROW .OR. MYCOL.GE.NPCOL ) GO TO 20
*
* Вычисление размеров локальных частей матрицы  A  на всех процессах
*
      NP = NUMROC( N, NB, MYROW, 0, NPROW )
      NQ = NUMROC( N, NB, MYCOL, 0, NPCOL )
*
* Инициализация дескрипторов для матриц A и Z
*
      CALL DESCINIT( DESCA, N, N, NB, NB, 0, 0, CTXT, MAX(1,NP), INFO)
      CALL DESCINIT( DESCZ, N, N, NB, NB, 0, 0, CTXT, MAX(1,NP), INFO)
*
* Распределение памяти  MEM
*
* Указатель на первый элемент массива собственных значений
      IPW = 1
* Указатель на первый элемент локальной части матрицы  A
      IPA = IPW + N
* Указатель на первый элемент массива собственных векторов
      IPZ = IPA + DESCA(9)*NQ
* Указатель на рабочий массив
      IPWORK = IPZ + DESCZ(9)*NQ
*
* Чтение из файла и распределение матрицы A
*
      CALL PDLAREAD('tdsyev25.dat', MEM(IPA), DESCA, 0, 0, MEM(IPWORK))
*
* Печать входных данных
*
      IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
         WRITE( NPR, FMT = 92 )
         WRITE( NPR, FMT = 98 )N, N, NB, NB
         WRITE( NPR, FMT = 97 )NPROW*NPCOL, NPROW, NPCOL
         WRITE( NPR, FMT = * ) ' Исходные данные '
         WRITE( NPR, FMT = 88 ) UPLO
         WRITE( NPR, FMT = * )'   IPW,   IPA,   IPZ,   IPWORK'
         WRITE( NPR, FMT = 94 )IPW, IPA, IPZ, IPWORK
         WRITE( NPR, FMT = * ) '  DESCA'
         WRITE( NPR, FMT = 85 )DESCA
         WRITE( NPR, FMT = * ) '  DESCZ'
         WRITE( NPR, FMT = 85 )DESCZ
         WRITE( NPR, FMT = * ) ' Матрица A'
      END IF

 94 FORMAT( / 4I7 /)
 85 FORMAT( / 9I5 /)
*  84 FORMAT( / ' MEM(IPWORK) = ', D16.5 /)
 88 FORMAT( / ' UPLO = ', A /)

      CALL PDLAPRNT( N, N, MEM(IPA), IA, JA, DESCA, 0, 0, 'A',
     $                                NPR, WORK1 )
*
* Вычисление собственных значений и собственных векторов
*
      CALL PDSYEV2( UPLO, N, MEM(IPA), IA, JA, DESCA, MEM(IPW),
     $                              MEM(IPZ), IZ, JZ, DESCZ, MEM(IPWORK), LWORK, INFO)
*
* Печать результатов
*
      IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
         WRITE( NPR, FMT = * ) ' Результаты '
         WRITE( NPR, FMT = 96 )INFO
*         WRITE( NPR, FMT = 84 )MEM(IPWORK)
         WRITE( NPR, FMT = * ) ' собственные числа'
*
         DO 10  I = 1, N
            PRINT *, ' W(', I, ')=', MEM( IPW+I-1 ), ';'
 10    CONTINUE
         WRITE( NPR, FMT = * ) ' собственные векторы'
      END IF
*
* Печать собственных векторов
*
      CALL PDLAPRNT( N, N, MEM(IPZ), IZ, JZ, DESCZ, 0, 0, 'Z', NPR,
     $                                WORK1 )
*
* Запись в файл результатов работы подпрограммы PDSYEV2.f
*
      IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
         OPEN( NOUT, FILE = 'tdsyev25.res', STATUS = 'UNKNOWN' )
         WRITE( NOUT, FMT = * ) ' собственные числа'
         DO 11  I = 1, N
            WRITE( NOUT, FMT = * ) MEM( IPW+I-1 )
 11    CONTINUE
         WRITE( NOUT, FMT = * ) ' собственные векторы'
      END IF
*
      CALL PDLAWRITE( 'tdsyev25.res', N, N, MEM(IPZ), 1, 1, DESCZ,
     $                                  0, 0, MEM(IPWORK) )
*
      CALL BLACS_GRIDEXIT( CTXT )
 20 CONTINUE

      CALL BLACS_EXIT( 0 )

 99 FORMAT( 'W=diag([', 4D16.12, //']);' )

 92 FORMAT( / 'Тест к подпрограмме PDSYEV2' )
 98 FORMAT( / ' Вычисление всех собственных значений и собственных',
     $  ' векторов матрицы  A,' //
     $      ' где A - симметричная матрица ',
     $      I2, ' на ', I2,', разбитая на блоки', I2, ' на ', I2, /)
 97 FORMAT( / ' Запуск на ', I2, ' процессах,',
     $       ' образующих решетку размером ', I2, ' на ', I2 /)
 96 FORMAT( / ' Значение INFO  = ', I6 /)
      STOP
      END

Файл с исходной матрицей:

           7           7
     0.10000D+01
     0.00000D+00
     0.00000D+00
     0.00000D+00
     0.00000D+00
     0.00000D+00
     0.10000D+01
     0.00000D+00
     0.10000D+01
     0.00000D+00
     0.00000D+00
     0.00000D+00
     0.00000D+00
     0.20000D+01
     0.00000D+00
     0.00000D+00
     0.10000D+01
     0.00000D+00
     0.00000D+00
     0.00000D+00
     0.30000D+01
     0.00000D+00
     0.00000D+00
     0.00000D+00
     0.10000D+01
     0.00000D+00
     0.00000D+00
     0.40000D+01
     0.00000D+00
     0.00000D+00
     0.00000D+00
     0.00000D+00
     0.10000D+01
     0.00000D+00
     0.50000D+01
     0.00000D+00
     0.00000D+00
     0.00000D+00
     0.00000D+00
     0.00000D+00
     0.10000D+01
     0.60000D+01
     0.10000D+01
     0.20000D+01
     0.30000D+01
     0.40000D+01
     0.50000D+01
     0.60000D+01
     0.70000D+01


Ниже приводится содержимое файла с результатами работы подпрограммы
PDSYEV2.

  собственные числа
  -6.00000000000000     
  0.999999999999999     
   1.00000000000000     
   1.00000000000000     
   1.00000000000000     
   1.00000000000000     
   14.0000000000000     
  собственные векторы
           7           7
     -0.845154254728516519D-01
     -0.169030850945703026D+00
     -0.253546276418554983D+00
     -0.338061701891406607D+00
     -0.422577127364258232D+00
     -0.507092552837109856D+00
      0.591607978309961480D+00
     -0.111381446794637101D+00
     -0.535086269928834346D+00
      0.819957700161629388D+00
     -0.330304485661413322D-01
     -0.155013470859046570D+00
     -0.618549942121308996D-01
      0.111022302462515654D-15
     -0.984609388674179908D+00
      0.126010706987393587D-01
     -0.987146606522586040D-01
      0.607396657301895460D-01
      0.921782793543347578D-01
      0.919501949235075333D-01
     -0.693889390390722838D-17
      0.000000000000000000D+00
     -0.179127814584487122D-01
      0.896727512209512106D-01
     -0.302109418637490967D+00
      0.802546495342590593D+00
     -0.506247915484824174D+00
      0.000000000000000000D+00
      0.000000000000000000D+00
     -0.710569201647523116D-02
      0.665136553040530758D-02
     -0.841076733874915994D+00
      0.233226858955070566D-01
      0.540325132243977158D+00
      0.000000000000000000D+00
     -0.845547991280188316D-01
      0.818044310655747453D+00
      0.459280136563631980D+00
     -0.143928223622574486D+00
     -0.219371577644360716D+00
     -0.209467241527044962D+00
     -0.222044604925031308D-15
     -0.620173672946042129D-01
     -0.124034734589208551D+00
     -0.186052101883812687D+00
     -0.248069469178416879D+00
     -0.310086836473021044D+00
     -0.372104203767625319D+00
     -0.806225774829854913D+00

Приложение 3_4

Пример составления головной программы
для целевой подпрограммы сингулярного разложения квадратной матрицы.

Требуется вычислить все сингулярные числа квадратной матрицы  A, которая имеет вид 

        |   1   2   3   4   5   |
        |   2   1   3   4   5   |
        |   3   3   1   4   5   |
        |   4   4   4   1   5   |
        |   5   5   5   5   1   |

Порядок матрицы N = 5.

Пусть матрица разбивается на блоки с размером блоков NB = 2.

Осуществляется пропуск программы на 4 процессах, которые образуют решетку 2 на 2.

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

      PROGRAM TDGESVD1
*
      include 'mpif.h'
* Именованные константы - параметры задачи
      INTEGER                        LWORK, MAXN, LDA, MAXPROCS, NOUT
      DOUBLE PRECISION   ZERO, MONE
      PARAMETER                 ( LDA = 10, LWORK = 50, MAXN = 100,
     $                                      MAXPROCS = 512, NOUT =6, ZERO = 0.0D+0,
     $                                      MONE = -1.0D+0 )
* Локальные переменные
*
      INTEGER            CTXT, I, IAM, INFO, MYCOL, MYROW, N, NB,
     $                           NPCOL, NPROCS, NPROW, IA, JA
*
* Локальные массивы
      INTEGER                        DESCA( 9 )
      DOUBLE PRECISION   A( LDA, LDA ), S( 5 ),
     $                                       WORK( LWORK ), WORK1( 5 )
*
* Внешние подпрограммы
      EXTERNAL     BLACS_EXIT, BLACS_GET, BLACS_GRIDEXIT,
     $                        BLACS_GRIDINFO, BLACS_GRIDINIT, BLACS_PINFO,
     $                        BLACS_SETUP, DESCINIT, PDMATINIT, PDLAPRNT,
     $                        PDGESVD1
*
* Выполняемые операторы
*
* Постановка задачи
*
      N = 5
      NB = 2
      NPROW = 2
      NPCOL = 2
      IA = 1
      JA = 1
*
* Инициализация пакета BLACS
*
      CALL  BLACS_PINFO( IAM, NPROCS )
      IF( ( NPROCS.LT.1 ) ) THEN
         CALL  BLACS_SETUP( IAM, NPROW*NPCOL )
      END IF
*
* Инициализация контекста BLACS'а и решетки процессов
*
      CALL  BLACS_GET( -1, 0, CTXT )
      CALL  BLACS_GRIDINIT( CTXT, 'R', NPROW, NPCOL )
      CALL  BLACS_GRIDINFO( CTXT, NPROW, NPCOL, MYROW, MYCOL )
*
* Если процесс не является частью данного контекста,
* переход на конец программы
*
      IF( MYROW.EQ.-1 ) GO TO 20
*
* Инициализация дескриптора для матрицы A
*
      CALL  DESCINIT( DESCA, N, N, NB, NB, 0, 0, CTXT, LDA, INFO )
*
* Построение квадратной матрицы A
*
      CALL  PDMATINIT( N, A, IA, JA, DESCA, INFO )
*
* Печать входных данных
*
      IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
         WRITE( NOUT, FMT = 92 )
         WRITE( NOUT, FMT = 98 )N, N, NB, NB
         WRITE( NOUT, FMT = 97 )NPROW*NPCOL, NPROW, NPCOL
         WRITE( NOUT, FMT = * ) ' Исходные данные '
         WRITE( NOUT, FMT = * ) '  DESCA'
         WRITE( NOUT, FMT = 85 )DESCA
         WRITE( NOUT, FMT = * ) ' Матрица  A'
      END IF
 85 FORMAT( / 9I5 /)

      CALL  PDLAPRNT( N, N, A, IA, JA, DESCA, 0, 0, 'A', NOUT, WORK1 )
*
* Вычисление всех сингулярных чисел квадратной матрицы  A
*
      CALL  PDGESVD1( N, A, IA, JA, DESCA, S, WORK, LWORK, INFO )
*
* Печать результатов
*
      IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
         WRITE( NOUT, FMT = * ) ' Результаты '
         WRITE( NOUT, FMT = 96 )INFO
         WRITE( NOUT, FMT = * ) ' Сингулярные числа '
         DO 10  I = 1, N
            PRINT *, ' S(', I, ')=', S( I ), ';'
 10    CONTINUE
      END IF
*
      CALL  BLACS_GRIDEXIT( CTXT )
 20 CONTINUE
      CALL  BLACS_EXIT( 0 )

 92 FORMAT( / 'Тест к подпрограмме PDGESVD1' )
 98 FORMAT( / ' Вычисление всех сингулярных чисел матрицы A,' //
     $      ' где A - квадратная матрица ',
     $      I2, ' на ', I2,',' /
     $       '      разбитая на блоки', I2, ' на ', I2, /)
 97 FORMAT( / ' Запуск на ', I2, ' процессах,',
     $       ' образующих решетку размером ', I2, ' на ', I2 /)
 96 FORMAT( / ' Значение INFO  = ', I6 /)
      STOP
      END
*
      SUBROUTINE  PDMATINIT( N, A, IA, JA, DESCA, INFO )
*
* PDMATINIT генерирует и распределяет матрицу A по решетке процессов
*
* Параметры
      INTEGER                        IA, INFO, JA, N
      INTEGER                        DESCA( * )
      DOUBLE PRECISION   A( * )
*
* Локальные переменные
      INTEGER                        I, J
      DOUBLE PRECISION   AIJ
*
* Внешние подпрограммы
      EXTERNAL            PDELSET
* Внутренние функции
      INTRINSIC          DBLE
*
* Выполняемые операторы
*
      INFO = 0
*
      IF( IA.NE.1 ) THEN
         INFO = -3
      ELSE IF( JA.NE.1 ) THEN
         INFO = -4
      END IF
*
      DO 2  J = 1, N
      DO 1  I = 1, N
      IF( I .EQ. J ) THEN
      CALL  PDELSET( A, I, J, DESCA, 1.0D0 )
      ELSE IF( I .GT. J ) THEN
      CALL  PDELSET( A, I, J, DESCA, DBLE(I) )
      ELSE IF( I .LT. J ) THEN
      CALL  PDELSET( A, I, J, DESCA, DBLE(J) )
      END IF
   1 CONTINUE
   2 CONTINUE
      RETURN
      END

  Результаты

  значение INFO  =   0

  сингулярные числа:

  S( 1 )=   17.2323289653599
  S( 2 )=   5.42527830414530
  S( 3 )=   3.51682585287666
  S( 4 )=   2.29022480833794
  S( 5 )=   0.999999999999999