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

Назначение

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

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

Подпрограмма PDGESVD3 выполняет сингулярное разложение (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  PDGESVD3 ( SV, N, A, IA, JA, DESCA, S,
                                             U, IU, JU, DESCU, VT, IVT, JVT, DESCVT,
                                             WORK, LWORK, INFO)

   Параметры 

   Версии 

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

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

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

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

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

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

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

      PROGRAM TDGESVD3
      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 )
*
      CHARACTER          SV
      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,
     $                              PDGESVD3
*
      N = 5
      NB = 2
      NPROW = 2
      NPCOL = 2
      IA = 1
      JA = 1
      IU = 1
      JU = 1
      IVT = 1
      JVT = 1
      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, 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  PDGESVD3( SV, 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;

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

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