## MB02SD

### LU factorization of an upper Hessenberg matrix

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

Purpose

```  To compute an LU factorization of an n-by-n upper Hessenberg
matrix H using partial pivoting with row interchanges.

```
Specification
```      SUBROUTINE MB02SD( N, H, LDH, IPIV, INFO )
C     .. Scalar Arguments ..
INTEGER           INFO, LDH, N
C     .. Array Arguments ..
INTEGER           IPIV(*)
DOUBLE PRECISION  H(LDH,*)

```
Arguments

Input/Output Parameters

```  N       (input) INTEGER
The order of the matrix H.  N >= 0.

H       (input/output) DOUBLE PRECISION array, dimension (LDH,N)
On entry, the n-by-n upper Hessenberg matrix to be
factored.
On exit, the factors L and U from the factorization
H = P*L*U; the unit diagonal elements of L are not stored,
and L is lower bidiagonal.

LDH     INTEGER
The leading dimension of the array H.  LDH >= max(1,N).

IPIV    (output) INTEGER array, dimension (N)
The pivot indices; for 1 <= i <= N, row i of the matrix
was interchanged with row IPIV(i).

```
Error Indicator
```  INFO    INTEGER
= 0:  successful exit;
< 0:  if INFO = -i, the i-th argument had an illegal
value;
> 0:  if INFO = i, U(i,i) is exactly zero. The
factorization has been completed, but the factor U
is exactly singular, and division by zero will occur
if it is used to solve a system of equations.

```
Method
```  The factorization has the form
H = P * L * U
where P is a permutation matrix, L is lower triangular with unit
diagonal elements (and one nonzero subdiagonal), and U is upper
triangular.

This is the right-looking Level 1 BLAS version of the algorithm

```
References
```  -

```
Numerical Aspects
```                             2
The algorithm requires 0( N ) operations.

```
```  None
```
Example

Program Text

```*     MB02SD 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          NMAX, NRHMAX
PARAMETER        ( NMAX = 20, NRHMAX = 20 )
INTEGER          LDB, LDH
PARAMETER        ( LDB = NMAX, LDH = NMAX )
INTEGER          LDWORK
PARAMETER        ( LDWORK = 3*NMAX )
INTEGER          LIWORK
PARAMETER        ( LIWORK = NMAX )
*     .. Local Scalars ..
DOUBLE PRECISION HNORM, RCOND
INTEGER          I, INFO, INFO1, J, N, NRHS
CHARACTER*1      NORM, TRANS
*     .. Local Arrays ..
DOUBLE PRECISION H(LDH,NMAX), B(LDB,NRHMAX), DWORK(LDWORK)
INTEGER          IPIV(NMAX), IWORK(LIWORK)
*     .. External Functions ..
DOUBLE PRECISION DLAMCH, DLANHS
EXTERNAL         DLAMCH, DLANHS
*     .. External Subroutines ..
EXTERNAL         DLASET, MB02RD, MB02SD, MB02TD
*     .. Executable Statements ..
*
WRITE ( NOUT, FMT = 99999 )
*     Skip the heading in the data file and read in the data.
READ ( NIN, FMT = '()' )
READ ( NIN, FMT = * ) N, NRHS, NORM, TRANS
IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
WRITE ( NOUT, FMT = 99994 ) N
ELSE
READ ( NIN, FMT = * ) ( ( H(I,J), J = 1,N ), I = 1,N )
IF ( NRHS.LT.0 .OR. NRHS.GT.NRHMAX ) THEN
WRITE ( NOUT, FMT = 99993 ) NRHS
ELSE
READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,NRHS ), I = 1,N )
IF ( N.GT.2 )
\$         CALL DLASET( 'Lower', N-2, N-2, ZERO, ZERO, H(3,1), LDH )
*           Compute the LU factorization of the upper Hessenberg matrix.
CALL MB02SD( N, H, LDH, IPIV, INFO )
*           Estimate the reciprocal condition number of the matrix.
HNORM = DLANHS( NORM, N, H, LDH, DWORK )
CALL MB02TD( NORM, N, HNORM, H, LDH, IPIV, RCOND, IWORK,
\$                   DWORK, INFO1 )
IF ( INFO.EQ.0 .AND. RCOND.GT.DLAMCH( 'Epsilon' ) ) THEN
*              Solve the linear system.
CALL MB02RD( TRANS, N, NRHS, H, LDH, IPIV, B, LDB, INFO )
*
WRITE ( NOUT, FMT = 99997 )
ELSE
WRITE ( NOUT, FMT = 99998 ) INFO
END IF
DO 10 I = 1, N
WRITE ( NOUT, FMT = 99996 ) ( B(I,J), J = 1,NRHS )
10          CONTINUE
WRITE ( NOUT, FMT = 99995 ) RCOND
END IF
END IF
STOP
*
99999 FORMAT (' MB02SD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from MB02SD = ',I2)
99997 FORMAT (' The solution matrix is ')
99996 FORMAT (20(1X,F8.4))
99995 FORMAT (/' Reciprocal condition number = ',D12.4)
99994 FORMAT (/' N is out of range.',/' N = ',I5)
99993 FORMAT (/' NRHS is out of range.',/' NRHS = ',I5)
END
```
Program Data
``` MB02SD EXAMPLE PROGRAM DATA
5    4      O      N
1.    2.    6.    3.    5.
-2.   -1.   -1.    0.   -2.
0.    3.    1.    5.    1.
0.    0.    2.    0.   -4.
0.    0.    0.    1.    4.
5.    5.    1.    5.
-2.    1.    3.    1.
0.    0.    4.    5.
2.    1.    1.    3.
-1.    3.    3.    1.
```
Program Results
``` MB02SD EXAMPLE PROGRAM RESULTS

The solution matrix is
0.0435   1.2029   1.6377   1.1014
1.0870  -4.4275  -5.5580  -2.9638
0.9130   0.7609  -0.1087   0.6304
-0.8261   2.4783   4.2174   2.7391
-0.0435   0.1304  -0.3043  -0.4348

Reciprocal condition number =   0.1554D-01
```