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

Назначение

Решение систем линейных алгебраических уравнений с прямоугольными матрицами полного ранга

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

Подпрограмма решает переопределенную или недоопределенную систему линейных уравнений с вещественной матрицей размера M на N  sub (A) = A ( IA : IA + M - 1, JA : JA + N - 1) или с ее транспонированной, с использованием QR - или LQ - разложения sub (A). Предполагается, что sub (A) имеет полный ранг.

Подпрограмма обеспечивает следующие возможности.

Если первый параметр TRANS = ' N ' и  m ≥ n (или  если TRANS = ' T ' и  m < n), то ищется решение переопределенной системы, т.е. решается задача наименьших квадратов

                      минимизировать  || sub(B) - sub(A) * X ||
            ( или минимизировать  || sub(B) - sub(A)T * X || )

Если TRANS = ' N ' и  m < n ( или  если TRANS = ' T ' и  m ≥ n), то ищется решение недоопределенной системы с минимальной нормой

                       sub(A) * X   =   sub(B)
             ( или sub(A)T * X  =   sub(B) )

  При этом
  sub(B) означает B( IB : IB + M - 1, JB : JB + NRHS - 1)     при TRANS = ' N '   и

  sub(B) означает B( IB : IB + N - 1, JB : JB + NRHS - 1)     при TRANS = ' T '.

Несколько векторов правых частей  b и векторов решений x могут быть обработаны за один вызов подпрограммы. При TRANS = ' N ' векторы решений хранятся в виде столбцов матрицы правых частей sub (B) размером N на NRHS, а при TRANS = ' T ' - в виде столбцов матрицы правых частей sub (B) размером M на NRHS.

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

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

         CALL  PDGELS ( TRANS, M, N, NRHS, A, IA, JA, DESCA, B, IB,
                                       JB, DESCB, WORK, LWORK, INFO)

   Параметры
TRANS - если TRANS = ' N ', решается система с матрицей sub (A),
если TRANS = ' T ', решается система с матрицей sub (A)T
(глобальный входной параметр, тип символьный);
M - число строк распределенной подматрицы sub (A); M ≥ 0 (глобальный входной параметр, тип целый);
N - число столбцов распределенной подматрицы sub (A); N ≥ 0 (глобальный входной параметр, тип целый);
NRHS - число правых частей, т.е. число столбцов распределенных подматриц sub (B) и X; NRHS ≥ 0 (глобальный входной параметр, тип целый);
A - указатель на локальную память, занимаемую массивом двойной точности с локальной размерностью ( LLD_A, LOCc ( JA + N - 1));
на входе - содержит матрицу A размера M на N;
если M ≥ N, на место sub (A) записывается информация о QR - разложении sub (A), выдаваемом подпрограммой PDGEQRF;
если M < N, на место sub (A) записывается информация о LQ - разложении sub (A), выдаваемом подпрограммой PDGELQF (локальный входной и локальный выходной параметр);
IA - номер строки в глобальном массиве A, который указывает на первую строку подматрицы sub (A) (глобальный входной параметр);
JA - номер столбца в глобальном массиве A, который указывает на первый столбец подматрицы sub (A) (глобальный входной параметр, тип целый);
DESCA - дескриптор распределенной матрицы A (одномерный массив длины DLEN_);
DESCA содержит информацию о размещении A в памяти; полное описание DESCA см. в разделе документации "Дескрипторы глобальных массивов" (глобальный и локальный входной параметр, тип целый);
B - указатель на локальную память, занимаемую массивом двойной точности с локальной размерностью ( LLD_B, LOCc ( JB + NRHS - 1));
на входе этот массив содержит локальные части распределенной матрицы B, которая состоит из векторов правых частей и которая хранится по столбцам;
sub (B) имеет размер M на NRHS, если TRANS = ' N ', и размер N на NRHS, если TRANS = ' T ';
на выходе на месте sub (B) записываются векторы решения, сохраняемые по столбцам:
если TRANS = ' N ' и M ≥ N, строки с 1 - ой по N - ю sub (B) содержат векторы решения задачи наименьших квадратов; остаточная сумма квадратов для решения в каждом столбце задается суммой квадратов элементов с номерами от (N + 1) - го до M - го в этом столбце;
если TRANS = ' N ' и M < N, строки с 1 - ой по N - ю sub (B) содержат векторы решений с минимальной нормой;
если TRANS = ' T ' и M ≥ N, строки с 1 - ой по M - ю sub (B) содержат векторы решений с минимальной нормой;
если TRANS = ' T ' и M < N, строки с 1 - ой по M - ю sub (B) содержат векторы решений задачи наименьших квадратов; остаточная сумма квадратов для решения в каждом столбце задается суммой квадратов элементов с номерами от (M + 1) - го до N - го в этом столбце
(локальный входной и локальный выходной параметр);
IB - номер строки в глобальном массиве B, указывающий первую строку подматрицы sub (B) (глобальный входной параметр, тип целый);
JB - номер столбца в глобальном массиве B, указывающий первый столбец sub (B) (глобальный входной параметр, тип целый);
DESCB - дескриптор распределенной матрицы B (одномерный массив длины DLEN_);
т.к. векторы решений X размещаются на месте векторов правых частей B, то в качестве 3-го и 4-го элементов массива DESCB следует указывать размеры max(M, N) и NRHS соответственно;
DESCB содержит информацию о размещении B в памяти; полное описание DESCB см. в разделе документации "Дескрипторы глобальных массивов" (глобальный и локальный входной параметр, тип целый);
WORK - одномерный рабочий массив двойной точности длины LWORK; на выходе - WORK (1) содержит минимальную требуемую величину LWORK (локальное рабочее пространство, локальный выходной параметр);
LWORK - размер рабочего пространства (массива) WORK;
минимальная требуемая величина LWORK может быть вычислена самой подпрограммой, если к ней обратиться со значением LWORK = - 1; при этом подпрограмма не производит никаких других вычислений, а требуемое значение LWORK возвращается в элементе массива WORK (1); при этом сообщений об ошибках не выдается (локальный или глобальный входной параметр);
INFO - целая переменная, диагностирующая результат работы подпрограммы (глобальный выходной параметр)
= 0 - успешное завершение работы;
< 0 - если i - ый фактический параметр является массивом и его j - ый элемент имеет недопустимое значение,
тогда INFO = - ( i * 100 + j ),
если i - ый фактический параметр является скаляром и имеет недопустимое значение, тогда INFO = - i 
   Версии
PSGELS -
PCGELS   
PZGELS     		Решение систем линейных алгебраических уравнений с прямоугольными матрицами полного ранга с использованием QR - или LQ - разложения для вещественных данных одинарной точности, комплексных данных одинарной точности и комплексных данных двойной точности соответственно 
   Вызываемые подпрограммы
PDGELQF - LQ - разложение прямоугольной матрицы методом отражений
PDGEQRF - QR - разложение прямоугольной матрицы методом отражений
PDORMLQ - умножение матрицы общего вида на ортогональную матрицу, полученную из LQ - разложения подпрограммой PDGELQF
PDORMQR - умножение матрицы общего вида на ортогональную матрицу, полученную из QR - разложения подпрограммой PDGEQRF
PDLABAD - вычисление квадратного корня из наименьшего и наибольшего представимого в машине числа, если модуль экспоненты выходит за допустимые пределы
PDLAMCH - вычисление машинных параметров для арифметики с плавающей запятой
PDLASET - внедиагональные элементы матрицы полагаются равными ALFA, а диагональные элементы - равными BETA
PDLASCL - умножение вещественной матрицы на вещественный скаляр
PDLANGE - вычисление значения 1 - ой нормы, нормы Фробениуса, бесконечной нормы или наибольшего абсолютного значения любого элемента треугольной матрицы общего вида 
   Замечания по использованию
  1. 

В подпрограммах PSGELS, PCGELS, PZGELS параметры A, B, WORK имеют тип REAL, COMPLEX и DOUBLE COMPLEX соответственно.

  2. 

Используются подпрограммы BLACS_GRIDINFO ( из пакета BLACS); LSAME, ILCM, INDXG2P, NUMROC, CHK1MAT, PCHK2MAT, PXERBLA ( из библиотеки ScaLAPACK_TOOLS);
PDTRSM ( из пакета PBLAS).

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

1. Решение переопределенной системы A X = B , где A - матрица размера 6 на 5.
Пусть матрица А системы имеет вид:

      36      -630        3360      -7560        7560
  -630    14700    -88200    211680   -220500
  3360  -88200    564480 -1411200   1512000
 -7560  211680 -1411200  3628800 -3969000
  7560 -220500  1512000 -3969000   4410000
 -2772    83160  -582120  1552320 -1746360

Вектор правых частей В имеет вид:

     463
 -13860
   97020
-258720
  291060
-116424

Матрица разбивается на блоки размеров 2 на 2.
Осуществляется пропуск программы на 6 процессах и используется двумерная решетка процессов 2 на 3.

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

      PROGRAM  TDGELS
      include 'mpif.h'
*
      INTEGER        DLEN_, IA, JA, IB, JB, N, NB, RSRC, CSRC,
     $                        MXLLDA, MXLLDB, NRHS, NBRHS, LWORK,
     $                        NOUT, MXLOCR, MXLOCC, MXRHSC
      PARAMETER    ( DLEN_ = 9, IA = 1, JA = 1, IB = 1, JB = 1,
     $                        M = 6, N = 5, MB = 2, NB = 2, RSRC = 0,
     $                        CSRC = 0, MXLLDA = 5, LWORK = 18,
*     $                        CSRC = 0, MXLLDA = 5, LWORK = -1,
     $                        MXLLDB = 5, NRHS = 1, NBRHS = 1, NOUT = 6,
     $                        MXLOCR = 5, MXLOCC = 4, MXRHSC = 1 )
*
      CHARACTER      TRANS
      PARAMETER      ( TRANS = 'N' )
*
      INTEGER            ICTXT, INFO, MYCOL, MYROW, NPCOL, NPROW
*
      INTEGER                        DESCA( DLEN_ ), DESCB( DLEN_ )
      DOUBLE PRECISION   A( MXLLDA, MXLOCC ), B( MXLLDB, MXRHSC ),
     $                                      WORK( LWORK ), WORK1( NB )
*     $                                    WORK( 18 ), WORK1( NB )
*
      EXTERNAL       BLACS_EXIT, BLACS_GRIDEXIT, BLACS_GRIDINFO,
     $                         DESCINIT, PDMATINIT, PDGELS,
     $                         PDLAPRNT, SL_INIT
*
      DATA               NPROW / 2 / , NPCOL / 3 /
*
      CALL  SL_INIT( ICTXT, NPROW, NPCOL )
      CALL  BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
*
      IF( MYROW.EQ.-1 ) GO TO 10
*
      CALL  DESCINIT( DESCA, M, N, MB, NB, RSRC, CSRC, ICTXT, MXLLDA,
     $               INFO )
      CALL  DESCINIT( DESCB, M, NRHS, NB, NBRHS, RSRC, CSRC, ICTXT,
     $               MXLLDB, INFO )
*
      CALL  PDMATINIT( M, N, A, IA, JA, DESCA, B, IB, JB, DESCB, INFO )
*
      CALL  PDGELS( TRANS, M, N, NRHS, A, IA, JA, DESCA,
     $              B, IB, JB, DESCB, WORK, LWORK, INFO )
*
      CALL  BLACS_GRIDEXIT( ICTXT )
 10 CONTINUE
      CALL  BLACS_EXIT( 0 )
      STOP
      END
*
      SUBROUTINE  PDMATINIT( M, N, A, IA, JA, DESCA, B, IB, JB, DESCB, INFO )
*
* PDMATINIT генерирует и распределяет матрицы A и B  по решетке процессов 2 x 3
*
      INTEGER            IA, IB, INFO, JA, JB, M, N
      INTEGER            DESCA( * ), DESCB( * )
      DOUBLE PRECISION   A( * ), B( * )
*
      INTEGER            I, J
      DOUBLE PRECISION   AIJ
*
      EXTERNAL           PDELSET
*
      INFO = 0
*
       CALL  PDELSET( A, 1, 1, DESCA, 36.0D0 )
       CALL  PDELSET( A, 2, 2, DESCA, 14700.0D0 )
       CALL  PDELSET( A, 3, 3, DESCA, 564480.0D0 )
       CALL  PDELSET( A, 4, 4, DESCA, 3628800.0D0 )
       CALL  PDELSET( A, 5, 5, DESCA, 4410000.0D0 )
*
      I = 6
      DO 3  J = 1, 5
         IF( J.EQ.1 ) AIJ  = -2772.0D0
         IF( J.EQ.2 ) AIJ  = 83160.0D0
         IF( J.EQ.3 ) AIJ  = -582120.0D0
         IF( J.EQ.4 ) AIJ  = 1552320.0D0
         IF( J.EQ.5 ) AIJ  = -1746360.0D0
         CALL  PDELSET( A, I, J, DESCA, AIJ )
   3 CONTINUE
*
      CALL  PDELSET( A, 1, 2 ,DESCA, -630.0D0 )
      CALL  PDELSET( A, 2, 1 ,DESCA, -630.0D0 )
      CALL  PDELSET( A, 1, 3 ,DESCA, 3360.0D0 )
      CALL  PDELSET( A, 3, 1 ,DESCA, 3360.0D0 )
      CALL  PDELSET( A, 1, 4 ,DESCA, -7560.0D0 )
      CALL  PDELSET( A, 4, 1 ,DESCA, -7560.0D0 )
      CALL  PDELSET( A, 1, 5 ,DESCA, 7560.0D0 )
      CALL  PDELSET( A, 5, 1 ,DESCA, 7560.0D0 )
      CALL  PDELSET( A, 2, 3 ,DESCA, -88200.0D0 )
      CALL  PDELSET( A, 3, 2 ,DESCA, -88200.0D0 )
      CALL  PDELSET( A, 2, 4 ,DESCA, 211680.0D0 )
      CALL  PDELSET( A, 4, 2 ,DESCA, 211680.0D0 )
      CALL  PDELSET( A, 2, 5 ,DESCA, -220500.0D0 )
      CALL  PDELSET( A, 5, 2 ,DESCA, -220500.0D0 )
      CALL  PDELSET( A, 3, 4 ,DESCA, -1411200.0D0 )
      CALL  PDELSET( A, 4, 3 ,DESCA, -1411200.0D0 )
      CALL  PDELSET( A, 3, 5 ,DESCA, 1512000.0D0 )
      CALL  PDELSET( A, 5, 3 ,DESCA, 1512000.0D0 )
      CALL  PDELSET( A, 4, 5 ,DESCA, -3969000.0D0 )
      CALL  PDELSET( A, 5, 4 ,DESCA, -3969000.0D0 )
*
      CALL  PDELSET( B, 1, 1, DESCB, 463.0D0 )
      CALL  PDELSET( B, 2, 1, DESCB, -13860.0D0 )
      CALL  PDELSET( B, 3, 1, DESCB, 97020.0D0 )
      CALL  PDELSET( B, 4, 1, DESCB, -258720.0D0 )
      CALL  PDELSET( B, 5, 1, DESCB, 291060.0D0 )
      CALL  PDELSET( B, 6, 1, DESCB, -116424.0D0 )
*
      RETURN
      END

Результаты:

  Значение   INFO  =   0

  Решение x

   x( 1, 1)  =  0.999999999944119811D+00
   x( 2, 1)  =  0.499999999981676879D+00
   x( 3, 1)  =  0.333333333325581516D+00
   x( 4, 1)  =  0.249999999996640215D+00
   x( 5, 1)  =  0.199999999998813599D+00

2. Решение недоопределенной системы A X = B , где A - матрица размера 5 на 6.
Пусть матрица А системы имеет вид:

      36      -630        3360      -7560        7560       -2772
  -630    14700    -88200    211680   -220500       83160
  3360  -88200    564480 -1411200   1512000    -582120
 -7560  211680 -1411200  3628800 -3969000    1552320
  7560 -220500  1512000 -3969000   4410000  -1746360

Вектор правых частей В имеет вид:

     53872.98
  -1542941.4
 10456462.8
-27223610.4
30065729.4

Матрица разбивается на блоки размеров 2 на 2.
Осуществляется пропуск программы на 6 процессах и используется двумерная решетка процессов 2 на 3.

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

      PROGRAM  TDGELS1
      include 'mpif.h'
*
      INTEGER        DLEN_, IA, JA, IB, JB, N, NB, RSRC, CSRC,
     $                        MXLLDA, MXLLDB, NRHS, NBRHS, LWORK,
     $                        NOUT, MXLOCR, MXLOCC, MXRHSC
      PARAMETER    ( DLEN_ = 9, IA = 1, JA = 1, IB = 1, JB = 1,
     $                        M = 5, N = 6, MB = 2, NB = 2, RSRC = 0,
     $                        CSRC = 0, MXLLDA = 5, LWORK = 23,
*     $                        CSRC = 0, MXLLDA = 5, LWORK = -1,
     $                        MXLLDB = 5, NRHS = 1, NBRHS = 1, NOUT = 6,
     $                        MXLOCR = 5, MXLOCC = 4, MXRHSC = 1 )
*
      CHARACTER      TRANS
      PARAMETER      ( TRANS = 'N' )
*
      INTEGER            ICTXT, INFO, MYCOL, MYROW, NPCOL, NPROW
*
      INTEGER                        DESCA( DLEN_ ), DESCB( DLEN_ )
      DOUBLE PRECISION   A( MXLLDA, MXLOCC ), B( MXLLDB, MXRHSC ),
     $                                      WORK( LWORK ), WORK1( NB )
*     $                                    WORK( 23 ), WORK1( NB )
*
      EXTERNAL       BLACS_EXIT, BLACS_GRIDEXIT, BLACS_GRIDINFO,
     $                         DESCINIT, PDMATINIT, PDGELS,
     $                         PDLAPRNT, SL_INIT
*
      DATA               NPROW / 2 / , NPCOL / 3 /
*
      CALL  SL_INIT( ICTXT, NPROW, NPCOL )
      CALL  BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
*
      IF( MYROW.EQ.-1 ) GO TO 10
*
      CALL  DESCINIT( DESCA, M, N, MB, NB, RSRC, CSRC, ICTXT, MXLLDA,
     $               INFO )
      CALL  DESCINIT( DESCB, N, NRHS, NB, NBRHS, RSRC, CSRC, ICTXT,
     $               MXLLDB, INFO )
*
      CALL  PDMATINIT( M, N, A, IA, JA, DESCA, B, IB, JB, DESCB, INFO )
*
      CALL  PDGELS( TRANS, M, N, NRHS, A, IA, JA, DESCA,
     $              B, IB, JB, DESCB, WORK, LWORK, INFO )
*
      CALL  BLACS_GRIDEXIT( ICTXT )
 10 CONTINUE
      CALL  BLACS_EXIT( 0 )
      STOP
      END
*
      SUBROUTINE  PDMATINIT( M, N, A, IA, JA, DESCA, B, IB, JB, DESCB, INFO )
*
* PDMATINIT генерирует и распределяет матрицы A и B  по решетке процессов 2 x 3
*
      INTEGER            IA, IB, INFO, JA, JB, M, N
      INTEGER            DESCA( * ), DESCB( * )
      DOUBLE PRECISION   A( * ), B( * )
*
      INTEGER            I, J
      DOUBLE PRECISION   AIJ
*
      EXTERNAL           PDELSET
*
      INFO = 0
*
       CALL  PDELSET( A, 1, 1, DESCA, 36.0D0 )
       CALL  PDELSET( A, 2, 2, DESCA, 14700.0D0 )
       CALL  PDELSET( A, 3, 3, DESCA, 564480.0D0 )
       CALL  PDELSET( A, 4, 4, DESCA, 3628800.0D0 )
       CALL  PDELSET( A, 5, 5, DESCA, 4410000.0D0 )
*
      J = 6
      DO 3  I = 1, 5
         IF( I.EQ.1 ) AIJ = -2772.0D0
         IF( I.EQ.2 ) AIJ = 83160.0D0
         IF( I.EQ.3 ) AIJ = -582120.0D0
         IF( I.EQ.4 ) AIJ = 1552320.0D0
         IF( I.EQ.5 ) AIJ = -1746360.0D0
         CALL  PDELSET( A, I, J, DESCA, AIJ )
   3 CONTINUE
*
      CALL  PDELSET( A, 1, 2 ,DESCA, -630.0D0 )
      CALL  PDELSET( A, 2, 1 ,DESCA, -630.0D0 )
      CALL  PDELSET( A, 1, 3 ,DESCA, 3360.0D0 )
      CALL  PDELSET( A, 3, 1 ,DESCA, 3360.0D0 )
      CALL  PDELSET( A, 1, 4 ,DESCA, -7560.0D0 )
      CALL  PDELSET( A, 4, 1 ,DESCA, -7560.0D0 )
      CALL  PDELSET( A, 1, 5 ,DESCA, 7560.0D0 )
      CALL  PDELSET( A, 5, 1 ,DESCA, 7560.0D0 )
      CALL  PDELSET( A, 2, 3 ,DESCA, -88200.0D0 )
      CALL  PDELSET( A, 3, 2 ,DESCA, -88200.0D0 )
      CALL  PDELSET( A, 2, 4 ,DESCA, 211680.0D0 )
      CALL  PDELSET( A, 4, 2 ,DESCA, 211680.0D0 )
      CALL  PDELSET( A, 2, 5 ,DESCA, -220500.0D0 )
      CALL  PDELSET( A, 5, 2 ,DESCA, -220500.0D0 )
      CALL  PDELSET( A, 3, 4 ,DESCA, -1411200.0D0 )
      CALL  PDELSET( A, 4, 3 ,DESCA, -1411200.0D0 )
      CALL  PDELSET( A, 3, 5 ,DESCA, 1512000.0D0 )
      CALL  PDELSET( A, 5, 3 ,DESCA, 1512000.0D0 )
      CALL  PDELSET( A, 4, 5 ,DESCA, -3969000.0D0 )
      CALL  PDELSET( A, 5, 4 ,DESCA, -3969000.0D0 )
*
      CALL  PDELSET( B, 1, 1, DESCB, 53872.98D0 )
      CALL  PDELSET( B, 2, 1, DESCB, -1542941.4D0 )
      CALL  PDELSET( B, 3, 1, DESCB, 10456462.8D0 )
      CALL  PDELSET( B, 4, 1, DESCB, -27223610.4D0 )
      CALL  PDELSET( B, 5, 1, DESCB, 30065729.4D0 )
*
      RETURN
      END

Результаты:

  Значение   INFO  =   0

  Решение x

 x( 1, 1) =  -0.202202646987226231D-10
 x( 2, 1) =  -0.700000000070026551D-01
 x( 3, 1) =   0.800000000003178613D+00
 x( 4, 1) =  -0.269999999999023110D+01
 x( 5, 1) =   0.350000000001390976D+01
 x( 6, 1) =  -0.153999999998352544D+01