MB04OD

QR factorization of a special structured block matrix (variant)

[Specification] [Arguments] [Method] [References] [Comments] [Example]

Purpose

  To calculate a QR factorization of the first block column and
  apply the orthogonal transformations (from the left) also to the
  second block column of a structured matrix, as follows
                       _   _
         [ R   B ]   [ R   B ]
    Q' * [       ] = [     _ ]
         [ A   C ]   [ 0   C ]
              _
  where R and R are upper triangular. The matrix A can be full or
  upper trapezoidal/triangular. The problem structure is exploited.

Specification
      SUBROUTINE MB04OD( UPLO, N, M, P, R, LDR, A, LDA, B, LDB, C, LDC,
     $                   TAU, DWORK )
C     .. Scalar Arguments ..
      CHARACTER         UPLO
      INTEGER           LDA, LDB, LDC, LDR, M, N, P
C     .. Array Arguments ..
      DOUBLE PRECISION  A(LDA,*), B(LDB,*), C(LDC,*), DWORK(*),
     $                  R(LDR,*), TAU(*)

Arguments

Mode Parameters

  UPLO    CHARACTER*1
          Indicates if the matrix A is or not triangular as follows:
          = 'U':  Matrix A is upper trapezoidal/triangular;
          = 'F':  Matrix A is full.

Input/Output Parameters
  N       (input) INTEGER                 _
          The order of the matrices R and R.  N >= 0.

  M       (input) INTEGER
          The number of columns of the matrices B and C.  M >= 0.

  P       (input) INTEGER
          The number of rows of the matrices A and C.  P >= 0.

  R       (input/output) DOUBLE PRECISION array, dimension (LDR,N)
          On entry, the leading N-by-N upper triangular part of this
          array must contain the upper triangular matrix R.
          On exit, the leading N-by-N upper triangular part of this
                                                     _
          array contains the upper triangular matrix R.
          The strict lower triangular part of this array is not
          referenced.

  LDR     INTEGER
          The leading dimension of array R.  LDR >= MAX(1,N).

  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
          On entry, if UPLO = 'F', the leading P-by-N part of this
          array must contain the matrix A. If UPLO = 'U', the
          leading MIN(P,N)-by-N part of this array must contain the
          upper trapezoidal (upper triangular if P >= N) matrix A,
          and the elements below the diagonal are not referenced.
          On exit, the leading P-by-N part (upper trapezoidal or
          triangular, if UPLO = 'U') of this array contains the
          trailing components (the vectors v, see Method) of the
          elementary reflectors used in the factorization.

  LDA     INTEGER
          The leading dimension of array A.  LDA >= MAX(1,P).

  B       (input/output) DOUBLE PRECISION array, dimension (LDB,M)
          On entry, the leading N-by-M part of this array must
          contain the matrix B.
          On exit, the leading N-by-M part of this array contains
                              _
          the computed matrix B.

  LDB     INTEGER
          The leading dimension of array B.  LDB >= MAX(1,N).

  C       (input/output) DOUBLE PRECISION array, dimension (LDC,M)
          On entry, the leading P-by-M part of this array must
          contain the matrix C.
          On exit, the leading P-by-M part of this array contains
                              _
          the computed matrix C.

  LDC     INTEGER
          The leading dimension of array C.  LDC >= MAX(1,P).

  TAU     (output) DOUBLE PRECISION array, dimension (N)
          The scalar factors of the elementary reflectors used.

Workspace
  DWORK   DOUBLE PRECISION array, dimension (MAX(N-1,M))

Method
  The routine uses N Householder transformations exploiting the zero
  pattern of the block matrix.  A Householder matrix has the form

                                  ( 1 )
     H  = I - tau *u *u',    u  = ( v ),
      i          i  i  i      i   (  i)

  where v  is a P-vector, if UPLO = 'F', or a min(i,P)-vector, if
         i
  UPLO = 'U'.  The components of v  are stored in the i-th column
                                  i
  of A, and tau  is stored in TAU(i).
               i
  In-line code for applying Householder transformations is used
  whenever possible (see MB04OY routine).

Numerical Aspects
  The algorithm is backward stable.

Further Comments
  None
Example

Program Text

*     MB04OD EXAMPLE PROGRAM TEXT.
*     Copyright (c) 2002-2017 NICONET e.V.
*
*     .. Parameters ..
      DOUBLE PRECISION ZERO
      PARAMETER        (ZERO  = 0.0D0 )
      INTEGER          NIN, NOUT
      PARAMETER        ( NIN = 5, NOUT = 6 )
      INTEGER          MMAX, NMAX, PMAX
      PARAMETER        ( MMAX = 20, NMAX = 20, PMAX = 20 )
      INTEGER          LDA, LDB, LDC, LDR
      PARAMETER        ( LDA = PMAX, LDB = NMAX, LDC = PMAX,
     $                   LDR = NMAX )
      INTEGER          LDWORK
      PARAMETER        ( LDWORK = MAX( NMAX-1,MMAX ) )
*     .. Local Scalars ..
      CHARACTER*1      UPLO
      INTEGER          I, J, M, N, P
*     .. Local Arrays ..
      DOUBLE PRECISION A(LDA,NMAX), B(LDB,MMAX), C(LDC,MMAX),
     $                 DWORK(LDWORK), R(LDR,NMAX), TAU(NMAX)
*     .. External Subroutines ..
      EXTERNAL         MB04OD
*     .. Intrinsic Functions ..
      INTRINSIC        MAX
*     .. Executable Statements ..
*
      WRITE ( NOUT, FMT = 99999 )
*     Skip the heading in the data file and read the data.
      READ ( NIN, FMT = '()' )
      READ ( NIN, FMT = * ) N, M, P, UPLO
      IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
         WRITE ( NOUT, FMT = 99993 ) N
      ELSE
         IF ( M.LT.0 .OR. M.GT.MMAX ) THEN
            WRITE ( NOUT, FMT = 99992 ) M
         ELSE
            IF ( P.LT.0 .OR. P.GT.PMAX ) THEN
               WRITE ( NOUT, FMT = 99991 ) P
            ELSE
               READ ( NIN, FMT = * ) ( ( R(I,J), J = 1,N ), I = 1,N )
               READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,P )
               READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,M ), I = 1,N )
               READ ( NIN, FMT = * ) ( ( C(I,J), J = 1,M ), I = 1,P )
*              Compute and apply QR factorization.
               CALL MB04OD( UPLO, N, M, P, R, LDR, A, LDA, B, LDB, C,
     $                      LDC,  TAU, DWORK )
*
               WRITE ( NOUT, FMT = 99997 )
               DO 40 I = 1, N
                  DO 20 J = 1, I-1
                     R(I,J) = ZERO
   20             CONTINUE
                  WRITE ( NOUT, FMT = 99996 ) ( R(I,J), J = 1,N )
   40          CONTINUE
               IF ( M.GT.0 ) THEN
                  WRITE ( NOUT, FMT = 99995 )
                  DO 60 I = 1, N
                     WRITE ( NOUT, FMT = 99996 ) ( B(I,J), J = 1,M )
   60             CONTINUE
                  IF ( P.GT.0 ) THEN
                     WRITE ( NOUT, FMT = 99994 )
                     DO 80 I = 1, P
                        WRITE ( NOUT, FMT = 99996 ) ( C(I,J), J = 1,M )
   80                CONTINUE
                  END IF
               END IF
            END IF
         END IF
      END IF
      STOP
*
99999 FORMAT (' MB04OD EXAMPLE PROGRAM RESULTS',/1X)
99997 FORMAT (' The updated matrix R is ')
99996 FORMAT (20(1X,F10.4))
99995 FORMAT (' The updated matrix B is ')
99994 FORMAT (' The updated matrix C is ')
99993 FORMAT (/' N is out of range.',/' N = ',I5)
99992 FORMAT (/' M is out of range.',/' M = ',I5)
99991 FORMAT (/' P is out of range.',/' P = ',I5)
      END
Program Data
 MB04OD EXAMPLE PROGRAM DATA
   3     2     2     F
   3.    2.    1.
   0.    2.    1.
   0.    0.    1.
   2.    3.    1.
   4.    6.    5.
   3.    2.
   1.    3.
   3.    2.
   1.    3.
   3.    2.
Program Results
 MB04OD EXAMPLE PROGRAM RESULTS

 The updated matrix R is 
    -5.3852    -6.6850    -4.6424
     0.0000    -2.8828    -2.0694
     0.0000     0.0000    -1.7793
 The updated matrix B is 
    -4.2710    -3.7139
    -0.1555    -2.1411
    -1.6021     0.9398
 The updated matrix C is 
     0.5850     1.0141
    -2.7974    -3.1162

Return to index