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

Назначение

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

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

Подпрограмма PDGESVD1 вычисляет сингулярные числа вещественной квадратной матрицы A

Матрица  A  приводится к двухдиагональной форме  A = U1 B V1T, где U1 и V1 являются ортогональными матрицами, а  B  является вещественной и верхней двухдиагональной.

Сингулярные числа матрицы A совпадают с сингулярными числами матрицы B.

Сингулярные числа матрицы B вычисляются с помощью  qd - алгоритма и помещаются в одномерном массиве двойной точности в убывающем порядке.

Матрица A должна быть задана в виде двумерного массива двойной точности. Все вычисления проводятся в режиме DOUBLE PRECISION.

Литература:
Fernando K.V. and Parlett B.N. Accurate singular values and differential qd algorithms// Numer.Math. 1994, Vol-67, No. 2, pp. 191-230.
Форсайт Дж., Малькольм М., Моулер К. Машинные методы математических вычислений. М.: Изд-во МИР, 1980.
http://www.netlib.org/scalapack/slug/index.html

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

         CALL  PDGESVD1 ( N, A, IA, JA, DESCA, S, WORK, LWORK, INFO)

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

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

PDGEBRD - приведение прямоугольной матрицы общего вида к вещественному двухдиагональному виду при помощи ортогональных преобразований;
PDLARED1D - перераспределяет одномерный массив, предполагая, что входной массив BYCOL является распределенным по строкам, и что все столбцы решетки процессов содержат ту же самую копию BYCOL; выходной массив BYALL будет идентичным на всех процессах и содержать весь массив;
PDLARED2D - перераспределяет одномерный массив в предположении, что входной массив BYROW является распределенным по столбцам, и что все строки решетки процессов содержат ту же самую копию BYROW; выходной массив BYALL будет идентичным на всех процессах и содержать весь массив;
PDLAMCH - вычисление машинных параметров для арифметики с плавающей запятой с удвоенной точностью;
PDLASCL - умножение вещественной матрицы на вещественный скаляр
DBDSQR - сингулярное разложение вещественной квадратной двухдиагональной (верхней или нижней) матрицы с использованием QR - алгоритма и вращения Ривеиса 
   Замечания по использованию
В подпрограмме PSGESVD1 параметры A, S имеют тип REAL.
Используются подпрограммы BLACS_GRIDINFO (из пакета BLACS), LSAME (из пакета BLAS), PXERBA (из пакета PBLAS), 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.

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

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

      PROGRAM TDGESVD1
      include 'mpif.h'
*
      INTEGER                   LWORK, MAXN, LDA, MAXPROCS, NOUT
      PARAMETER            ( LDA = 100, LWORK = 50, 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
*
      INTEGER                        DESCA( 9 )
      DOUBLE PRECISION   A( LDA, LDA ), S( 5 ),
     $                                       WORK( LWORK ), WORK1( 5 )
*     $                                       WORK( 50 ), 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
*
      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 )
*
* Построение матрицы A
*
      CALL  PDMATINIT( N, A, IA, JA, DESCA, INFO )
*
* Вычисление всех сингулярных чисел квадратной матрицы  A
*
      CALL  PDGESVD1( N, A, IA, JA, DESCA, S, 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;