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

Назначение

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

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

Подпрограмма PDGESVD5 выполняет сингулярное разложение (SVD - singular value decomposition) матрицы A размера M на N в виде 

                   A = U * D * VT  ,

  где
             D - прямоугольная матрица размера M на N, у которой
                   все элементы, кроме диагональных dk k,  k = 1, 2, ..., min(M,N), равны нулю;

        U, V - ортогональные матрицы порядка M и N соответственно.

Вычисляются левые и правые сингулярные векторы.

Диагональные элементы матрицы D являются сингулярными числами матрицы A, а столбцы матриц U и V являются соответствующими правыми и левыми сингулярными векторами соответственно.

Сингулярные числа помещаются в одномерном массиве двойной точности в убывающем порядке и вычисляются только первые min (M, N) столбцов матрицы U и строк матрицы VT = VT.
Матрицы A, U и VT должны быть заданы в виде двумерных массивов двойной точности. Все вычисления проводятся в режиме DOUBLE PRECISION.

Литература:
Форсайт Дж., Малькольм М., Моулер К. Машинные методы математических вычислений. М.: Изд-во МИР, 1980.
http://www.netlib.org/scalapack/slug/index.html

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

         CALL  PDGESVD5 ( M, N, A, IA, JA, DESCA, S,
                                            U, IU, JU, DESCU, VT, IVT, JVT, DESCVT,
                                            WORK, LWORK, INFO)

   Параметры
M - число строк входной матрицы A, (глобальный входной параметр, тип целый);
N - число столбцов входной матрицы A, (глобальный входной параметр, тип целый);
A - массив двойной точности глобальной размерности (M, N), локальной размерности (MP, NQ), распределенный блочно - циклическим образом (см.),
MP - число локальных строк в матрицах A и U;
NQ - число локальных столбцов в матрицах A и VT;
на выходе содержимое  A не сохраняется
(локальный входной параметр, локальное рабочее пространство);
IA - номер строки в глобальном массиве A, указывающий на первую строку подматрицы sub (A) (глобальный входной параметр, тип целый);
JA - номер столбца в глобальном массиве A, указывающий на первый столбец подматрицы sub (A) (глобальный входной параметр, тип целый);
DESCA - дескриптор распределенной матрицы A (одномерный массив длины DLEN_); подробное описание см. в разделе документации "Дескрипторы глобальных массивов" (глобальный и локальный входной параметр, тип целый);
S - одномерный массив двойной точности длины SIZE = min (M, N); содержит сингулярные числа матрицы A, расположенные по убыванию т.е. S ( i ) і S ( i + 1); (глобальный выходной параметр);
U - массив двойной точности глобальной размерности (M, SIZE), локальной размерности (MP, SIZEQ),
SIZEQ - число локальных столбцов в массиве U;
MP - число локальных строк в массивах  A и U;
на выходе в первых min (M, N) столбцах U содержатся правые сингулярные векторы (локальный выходной параметр);
IU - номер строки в глобальном массиве U, указывающий на первую строку подматрицы sub (U) (глобальный входной параметр, тип целый);
JU - номер столбца в глобальном массиве U, указывающий на первый столбец подматрицы sub (U) (глобальный входной параметр, тип целый);
DESCU - дескриптор распределенной матрицы U (одномерный массив длины DLEN_); подробное описание см.в разделе документации "Дескрипторы глобальных массивов" (глобальный и локальный входной параметр, тип целый);
VT - массив двойной точности глобальной размерности (SIZE, N), локальной размерности (SIZEP, NQ),
SIZEP - число локальных строк в массиве VT;
NQ - число локальных столбцов в массивах A и VT;
на выходе в первых SIZE строках VT содержатся левые сингулярные векторы (локальный выходной параметр);
IVT - номер строки в глобальном массиве VT, указывающий на первую строку подматрицы sub (VT) (глобальный входной параметр, тип целый);
JVT - номер столбца в глобальном массиве VT, указывающий на первый столбец подматрицы sub (VT) (глобальный входной параметр, тип целый);
DESCVT - дескриптор распределенной матрицы VT (одномерный массив длины DLEN_); подробное описание см.в разделе документации "Дескрипторы глобальных массивов" (глобальный и локальный входной параметр, тип целый);
WORK - одномерный рабочий массив двойной точности длины LWORK; на выходе (если INFO = 0) WORK (1) содержит минимальную требуемую величину LWORK (локальное рабочее пространство, локальный выходной параметр);
LWORK - размер рабочего пространства (массива) WORK; минимальная требуемая величина LWORK может быть вычислена самой подпрограммой, если к ней обратиться со значением LWORK = - 1; при этом подпрограмма не производит никаких других вычислений, а требуемое значение LWORK возвращается в элементе массива WORK (1); при этом, сообщений об ошибках не выдается (локальный входной параметр);
INFO - целая переменная, диагностирующая результат работы подпрограммы (глобальный выходной параметр):
= 0 - успешное завершение работы;
< 0 - если i-й фактический параметр является массивом и его j-й элемент имеет недопустимое значение, то INFO = -( i*100 + j ),
если i-й фактический параметр является скаляром и имеет недопустимое значение, то INFO = -i;
> 0 - если подпрограмма DBDSQR не обеспечивает сходимость;
если INFO = MIN(M, N)+1, то в этом случае PDGESVD5 не гарантирует точность полученных результатов. 
   Версии
PSGESVD5 - сингулярное разложение прямоугольной матрицы размера M на N с вычислением всех сингулярных чисел и сингулярных векторов для вещественных данных одинарной точности 
   Вызываемые подпрограммы

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

PDGEBRD - приведение прямоугольной матрицы общего вида к вещественному двухдиагональному виду при помощи ортогональных преобразований;
PDGEMR2D - перераспределение элементов матриц по процессам блочно-циклическим образом;
PDORMBR - умножение матрицы общего вида на ортогональную матрицу преобразований, полученную при приведении матрицы A к двухдиагональному виду подпрограммой PDGEBRD;
PDLARED1D - перераспределяет одномерный массив, предполагая, что входной массив BYCOL является распределенным по строкам, и что все столбцы решетки процессов содержат ту же самую копию BYCOL; выходной массив BYALL будет идентичным на всех процессах и содержать весь массив;
PDLARED2D - перераспределяет одномерный массив в предположении, что входной массив BYROW является распределенным по столбцам и что все строки решетки процессов содержат ту же самую копию BYROW; выходной массив BYALL будет идентичным на всех процессах и содержать весь массив;
PDLAMCH - вычисление машинных параметров для арифметики с плавающей запятой с удвоенной точностью;
PDLASET - присвоение всем внедиагональным элементам матрицы вещественного значения двойной точности равного  a, а всем диагональным элементам - вещественного значения двойной точности равного  b;
PDLASCL - умножение вещественной матрицы на вещественный скаляр
DBDSQR - сингулярное разложение вещественной квадратной двухдиагональной (верхней или нижней) матрицы с использованием QR - алгоритма и вращения Ривеиса 
   Замечания по использованию
В подпрограмме PSGESVВ5 параметры A, U, VT, S имеют тип REAL.
Используются подпрограммы: BLACS_GET, BLACS_GRIDEXIT, BLACS_GRIDINFO, BLACS_GRIDINIT (из пакета BLACS),
LSAME (из пакета BLAS), PXERBA (из пакета PBLAS), ICEIL, CHK1MAT, PCHK1MAT, NUMROC, (из библиотеки ScaLAPACK_TOOLS), PDLACPY (из пакета ScaLAPACK).

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

Вычисление всех сингулярных чисел и сингулярных векторов матрицы A, которая имеет вид 

        |   1    6   11   |
        |   2    7   12   |
        |   3    8   13   |
        |   4    9   14   |
        |   5  10   15   |

Размерность матрицы M = 5, N = 3.

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

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

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

Матрица A инициализируется посредством обращения к подпрограмме PDELSET из библиотеки ScaLAPACK_TOOLS. 

      PROGRAM TDGESVD5
      include 'mpif.h'
*
      INTEGER                   LWORK, MAXN, LDA, MAXPROCS, NOUT
      PARAMETER            ( LDA = 100, LWORK = 52, MAXN = 100,
*      PARAMETER          ( LDA = 100, LWORK = -1, MAXN = 100,
     $                                   MAXPROCS = 512, NOUT =6 )
*
      INTEGER            CTXT, I, IAM, INFO, MYCOL, MYROW, M, N, NB,
     $                           NPCOL, NPROCS, NPROW, IA, JA, JU, IU, IVT, JVT, SIZE
*
      INTEGER                       DESCA( 9 ), DESCU( 9 ), DESCVT( 9 )
      DOUBLE PRECISION   A( LDA, LDA ), S( 3 ), U(LDA,LDA), VT(LDA,LDA),
*                                           U(3,2), VT(2,2),
     $                                      WORK( LWORK ), WORK1( 5 )
*     $                                      WORK( 52 ), WORK1( 5 )
*
      EXTERNAL           BLACS_EXIT, BLACS_GET, BLACS_GRIDEXIT,
     $                              BLACS_GRIDINFO, BLACS_GRIDINIT, BLACS_PINFO,
     $                              BLACS_SETUP, DESCINIT, PDMATINIT, PDLAPRNT, PDGESVD5
*
      N = 3
      M = 5
      NB = 2
      NPROW = 2
      NPCOL = 2
      IA = 1
      JA = 1
      IU = 1
      JU = 1
      IVT = 1
      JVT = 1
      SIZE = MIN( M, N )
*
      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, M, N, NB, NB, 0, 0, CTXT, LDA, INFO )
      CALL  DESCINIT( DESCU, M, SIZE, NB, NB, 0, 0, CTXT, LDA, INFO )
      CALL  DESCINIT( DESCVT, SIZE, N, NB, NB, 0, 0, CTXT, LDA, INFO )
*
* Построение матрицы A
*
      CALL  PDMATINIT( M, N, A, IA, JA, DESCA, INFO )
*
* Вычисление всех сингулярных чисел и векторов матрицы  A
*
      CALL  PDGESVD5( M, N, A, IA, JA, DESCA, S, U, IU, JU, DESCU,
     $                                 VT, IVT, JVT, DESCVT, WORK, LWORK, INFO )
*
      CALL  BLACS_GRIDEXIT( CTXT )

 20 CONTINUE
      CALL  BLACS_EXIT( 0 )
      STOP
      END
*
      SUBROUTINE  PDMATINIT( M, N, A, IA, JA, DESCA, INFO )
*
* PDMATINIT генерирует и распределяет матрицу A по решетке процессов
*
* Параметры
      INTEGER                        IA, INFO, JA, M, 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, M
      AIJ = DBLE(I + (J-1)*M)
      CALL  PDELSET( A, I, J, DESCA, AIJ)
   1 CONTINUE
   2 CONTINUE
      RETURN
      END

Результаты:

 Значение INFO  =  0

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

  S( 1 ) = 35.12722;
  S( 2 ) =  2.46539;
  S( 3 ) =  3.1655 E-015;

 Сингулярные векторы:

  Матрица U

 U( 1, 1) =  -0.3545571D+00
 U( 2, 1) =  -0.3986964D+00
 U( 3, 1) =  -0.4428357D+00
 U( 4, 1) =  -0.4869750D+00
 U( 5, 1) =  -0.5311143D+00

 U( 1, 2) =  -0.6886866D+00
 U( 2, 2) =  -0.3755545D+00
 U( 3, 2) =  -0.6242242D-01
 U( 4, 2) =   0.2507097D+00
 U( 5, 2) =   0.5638418D+00

 U( 1, 3) =   0.5328844D+00
 U( 2, 3) =  -0.5695110D+00
 U( 3, 3) =   0.5727210D-03
 U( 4, 3) =  -0.4241502D+00
 U( 5, 3) =   0.4602041D+00

  Матрица VT

 VT( 1, 1) =  -0.2016649D+00
 VT( 2, 1) =   0.8903171D+00
 VT( 3, 1) =  -0.4082483D+00

 VT( 1, 2) =  -0.5168305D+00
 VT( 2, 2) =   0.2573316D+00
 VT( 3, 2) =   0.8164966D+00

 VT( 1, 3) =  -0.8319961D+00
 VT( 2, 3) =  -0.3756539D+00
 VT( 3, 3) =  -0.4082483D+00