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

Назначение

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

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

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

                   A = U * D * VT  ,

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

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

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

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

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

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

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

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

   Параметры
N - порядок исходной матрицы A, (глобальный входной параметр, тип целый);
A - массив двойной точности глобальной размерности (N, N), локальной размерности (MP, NQ), распределенный блочно - циклическим образом (см.),
MP - число локальных строк в матрицах A и U;
NQ - число локальных столбцов в матрицах A и VT;
на выходе содержимое  A не сохраняется
(локальный входной параметр, локальное рабочее пространство);
IA - номер строки в глобальном массиве A, указывающий на первую строку подматрицы sub (A) (глобальный входной параметр, тип целый);
JA - номер столбца в глобальном массиве A, указывающий на первый столбец подматрицы sub (A) (глобальный входной параметр, тип целый);
DESCA - дескриптор распределенной матрицы A (одномерный массив длины DLEN_); подробное описание см. в разделе документации "Дескрипторы глобальных массивов" (глобальный и локальный входной параметр, тип целый);
S - одномерный массив двойной точности длины N; содержит сингулярные числа матрицы A, расположенные по убыванию т.е. S ( i ) ≥ S ( i + 1); (глобальный выходной параметр);
U - массив двойной точности глобальной размерности (N, N), локальной размерности (MP, SIZEQ),
SIZEQ - число локальных столбцов в массиве U;
MP - число локальных строк в массивах  A и U;
на выходе в N столбцах U содержатся правые сингулярные векторы (локальный выходной параметр);
IU - номер строки в глобальном массиве U, указывающий на первую строку подматрицы sub (U) (глобальный входной параметр, тип целый);
JU - номер столбца в глобальном массиве U, указывающий на первый столбец подматрицы sub (U) (глобальный входной параметр, тип целый);
DESCU - дескриптор распределенной матрицы U (одномерный массив длины DLEN_); подробное описание см.в разделе документации "Дескрипторы глобальных массивов" (глобальный и локальный входной параметр, тип целый);
VT - массив двойной точности глобальной размерности (N, N), локальной размерности (SIZEP, NQ),
SIZEP - число локальных строк в массиве VT;
NQ - число локальных столбцов в массивах A и VT;
на выходе в N строках 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 = N+1, то в этом случае PDGESVD2 не гарантирует точность полученных результатов. 
   Версии
PSGESVD2 - сингулярное разложение квадратной матрицы порядка N с вычислением всех сингулярных чисел и сингулярных векторов для вещественных данных одинарной точности 
   Вызываемые подпрограммы

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

PDGEBRD - приведение прямоугольной матрицы общего вида к вещественному двухдиагональному виду при помощи ортогональных преобразований;
PDGEMR2D - перераспределение элементов матриц по процессам блочно-циклическим образом;
PDORMBR - умножение матрицы общего вида на ортогональную матрицу преобразований, полученную при приведении матрицы A к двухдиагональному виду подпрограммой PDGEBRD;
PDLARED1D - перераспределяет одномерный массив, предполагая, что входной массив BYCOL является распределенным по строкам, и что все столбцы решетки процессов содержат ту же самую копию BYCOL; выходной массив BYALL будет идентичным на всех процессах и содержать весь массив;
PDLARED2D - перераспределяет одномерный массив в предположении, что входной массив BYROW является распределенным по столбцам и что все строки решетки процессов содержат ту же самую копию BYROW; выходной массив BYALL будет идентичным на всех процессах и содержать весь массив;
PDLAMCH - вычисление машинных параметров для арифметики с плавающей запятой с удвоенной точностью;
PDLASET - присвоение всем внедиагональным элементам матрицы вещественного значения двойной точности равного  a, а всем диагональным элементам - вещественного значения двойной точности равного  b;
PDLASCL - умножение вещественной матрицы на вещественный скаляр
DBDSQR - сингулярное разложение вещественной квадратной двухдиагональной (верхней или нижней) матрицы с использованием QR - алгоритма и вращения Ривеиса 
   Замечания по использованию
В подпрограмме PSGESVD2 параметры 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   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.

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

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

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

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

  Матрица U

 U( 1, 1) =  -0.404652434305128983D+00
 U( 2, 1) =  -0.404652434305129094D+00
 U( 3, 1) =  -0.425692654396723935D+00
 U( 4, 1) =  -0.465693233173190169D+00
 U( 5, 1) =  -0.523859133156270529D+00

 U( 1, 2) =  -0.389791547190647403D+00
 U( 2, 2) =  -0.389791547190647347D+00
 U( 3, 2) =  -0.275993043075253353D+00
 U( 4, 2) =   0.454476316088957971D-01
 U( 5, 2) =   0.786058173694953566D+00

 U( 1, 3) =   0.247683855301384148D+00
 U( 2, 3) =   0.247683855301383815D+00
 U( 3, 3) =   0.843929575152923550D-01
 U( 4, 3) =  -0.874087946343403699D+00
 U( 5, 3) =   0.325812072251903173D+00

 U( 1, 4) =  -0.350673159755202413D+00
 U( 2, 4) =  -0.350673159755202524D+00
 U( 3, 4) =   0.857607971563088434D+00
 U( 4, 4) =  -0.130516617390331158D+00
 U( 5, 4) =  -0.391235191545428740D-01

 U( 1, 5) =  -0.707106781186547462D+00
 U( 2, 5) =   0.707106781186547906D+00
 U( 3, 5) =   0.138777878078144568D-15
 U( 4, 5) =  -0.222044604925031308D-15
 U( 5, 5) =  -0.555111512312578270D-16

  Матрица VT

 VT( 1, 1) =  -0.404652434305128816D+00
 VT( 2, 1) =   0.389791547190647181D+00
 VT( 3, 1) =  -0.247683855301384176D+00
 VT( 4, 1) =   0.350673159755202302D+00
 VT( 5, 1) =   0.707106781186547684D+00

 VT( 1, 2) =  -0.404652434305128872D+00
 VT( 2, 2) =   0.389791547190647458D+00
 VT( 3, 2) =  -0.247683855301383982D+00
 VT( 4, 2) =   0.350673159755202191D+00
 VT( 5, 2) =  -0.707106781186547684D+00

 VT( 1, 3) =  -0.425692654396723880D+00
 VT( 2, 3) =   0.275993043075253464D+00
 VT( 3, 3) =  -0.843929575152923134D-01
 VT( 4, 3) =  -0.857607971563088434D+00
 VT( 5, 3) =   0.555111512312578270D-16

 VT( 1, 4) =  -0.465693233173190002D+00
 VT( 2, 4) =  -0.454476316088956445D-01
 VT( 3, 4) =   0.874087946343404032D+00
 VT( 4, 4) =   0.130516617390330936D+00
 VT( 5, 4) =   0.000000000000000000D+00

 VT( 1, 5) =  -0.523859133156270307D+00
 VT( 2, 5) =  -0.786058173694953566D+00
 VT( 3, 5) =  -0.325812072251903229D+00
 VT( 4, 5) =   0.391235191545432626D-01
 VT( 5, 5) =  -0.499600361081320443D-15