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

Назначение

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

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

Подпрограмма PDGESVD6 выполняет сингулярное разложение (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 = V^T.
Матрицы A, U и VT должны быть заданы в виде двумерных массивов двойной точности. Все вычисления проводятся в режиме DOUBLE PRECISION.

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

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

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

   Параметры 

   Версии 

   Вызываемые подпрограммы 

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

   Замечания по использованию  

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

Вычисление сингулярных чисел и левых сингулярных векторов матрицы 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.

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

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

      PROGRAM TDGESVD6
      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 )
*
      CHARACTER          SV
      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,
     $                              PDGESVD6
*
      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 )
      SV = 'L'
*
      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  PDGESVD6( SV, 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 ;


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

  Матрица 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