◆ Zsytrf()

UDU^T or LDL^T factorization of a complex symmetric matrix

Purpose: This routine computes the factorization of a complex symmetric matrix A using the Bunch-Kaufman diagonal pivoting method. The form of the factorization is
A = U*D*U^T or A = L*D*L^T

where U (or L) is a product of permutation and unit upper (lower) triangular matrices, and D is symmetric and block diagonal with 1 x 1 and 2 x 2 diagonal blocks.
This is the blocked version of the algorithm, calling Level 3 BLAS.

Parameters

[in]	Uplo	= "U": Upper triangle of A is stored. = "L": Lower triangle of A is stored.
[in]	N	Order of the matrix A. (N >= 0) (If N = 0, returns without computation)
[in,out]	A()	Array A(LA1 - 1, LA2 - 1) (LA1 >= N, LA2 >= N) [in] N x N symmetric matrix A. The upper or lower triangular part is to be referenced in accordance with Uplo. [out] The block diagonal matrix D and the multipliers used to obtain the factor U or L (see below for further details).
[out]	IPiv()	Array IPiv(LIPiv - 1) (LIPiv >= N) Details of the interchanges and the block structure of D. If IPiv(k-1) > 0, then rows and columns k and IPiv(k-1) were interchanged, and k-th diagonal of D is a 1 x 1 diagonal block. If Uplo = "U" and IPiv(k-1) = IPiv(k-2) < 0, then rows and columns k-1 and -IPiv(k-1) were interchanged and (k-1)-th diagonal of D is a 2 x 2 diagonal block. If Uplo = "L" and IPiv(k-1) = IPiv(k) < 0, then rows and columns k+1 and -IPiv(k-1) were interchanged and k-th diagonal of D is a 2 x 2 diagonal block.
[out]	Info	= 0: Successful exit. = -1: The argument Uplo had an illegal value. (Uplo <> "U" nor "L") = -2: The argument N had an illegal value. (N < 0) = -3: The argument A() is invalid. = -4: The argument IPiv() is invalid. = i > 0: The i-th diagonal element of D is exactly zero. The factorization has been completed, but the block diagonal matrix D is exactly singular, and division by zero will occur if it is used to solve a system of equations.

Further Details: If Uplo = "U", then A = U*D*U^T, where
U = P(n)*U(n)* ... *P(k)U(k)* ...,

i.e., U is a product of terms P(k)*U(k), where k decreases from n to 1 in steps of 1 or 2, and D is a block diagonal matrix with 1 x 1 and 2 x 2 diagonal blocks D(k). P(k) is a permutation matrix as defined by IPiv(k-1), and U(k) is a unit upper triangular matrix, such that if the diagonal block D(k) is of order s (s = 1 or 2), then
( I v 0 ) k-s

U(k) = ( 0 I 0 ) s

( 0 0 I ) n-k

k-s s n-k

If s = 1, D(k) overwrites A(k-1, k-1), and v overwrites A(0 to k-2, k-1).
If s = 2, the upper triangle of D(k) overwrites A(k-2, k-2), A(k-2, k-1), and A(k-1, k-1), and v overwrites A(0 to k-3, k-2 to k-1).

If Uplo = "L", then A = L*D*L^T, where
L = P(1)*L(1)* ... *P(k)*L(k)* ...,

i.e., L is a product of terms P(k)*L(k), where k increases from 1 to n in steps of 1 or 2, and D is a block diagonal matrix with 1 x 1 and 2 x 2 diagonal blocks D(k). P(k) is a permutation matrix as defined by IPiv(k-1), and L(k) is a unit lower triangular matrix, such that if the diagonal block D(k) is of order s (s = 1 or 2), then
( I 0 0 ) k-1

L(k) = ( 0 I 0 ) s

( 0 v I ) n-k-s+1

k-1 s n-k-s+1

If s = 1, D(k) overwrites A(k-1, k-1), and v overwrites A(k to n-1, k-1).
If s = 2, the lower triangle of D(k) overwrites A(k-1, k-1), A(k-1, k), and A(k, k), and v overwrites A(k+1 to n-1, k-1 to k).

Example Program: Solve the system of linear equations Ax = B and estimate the reciprocal of the condition number (RCond) of A, where
( 0.20-0.11i -0.93-0.32i -0.80-0.92i )

A = ( -0.93-0.32i 0.81+0.37i -0.29+0.86i )

( -0.80-0.92i -0.29+0.86i 0.64+0.51i )

( 1.1120-1.0248i )

B = ( -1.5297-0.7781i )

( -0.4965-0.6057i )

Sub Ex_Zsytrf()

Const N As Long = 3

Dim A(N - 1, N - 1) As Complex, B(N - 1) As Complex, IPiv(N - 1) As Long

Dim ANorm As Double, RCond As Double, Info As Long

A(0, 0) = Cmplx(0.2, -0.11)

A(1, 0) = Cmplx(-0.93, -0.32): A(1, 1) = Cmplx(0.81, 0.37)

A(2, 0) = Cmplx(-0.8, -0.92): A(2, 1) = Cmplx(-0.29, 0.86): A(2, 2) = Cmplx(0.64, 0.51)

B(0) = Cmplx(1.112, -1.0248): B(1) = Cmplx(-1.5297, -0.7781): B(2) = Cmplx(-0.4965, -0.6057)

ANorm = Zlansy("1", "L", N, A())

Call Zsytrf("L", N, A(), IPiv(), Info)

If Info = 0 Then Call Zsytrs("L", N, A(), IPiv(), B(), Info)

If Info = 0 Then Call Zsycon("L", N, A(), IPiv(), ANorm, RCond, Info)

Debug.Print "X =",

Debug.Print Creal(B(0)), Cimag(B(0)), Creal(B(1)), Cimag(B(1)), Creal(B(2)), Cimag(B(2))

Debug.Print "RCond =", RCond

Debug.Print "Info =", Info

End Sub

Example Results: X = 0.71 0.59 -0.15 0.19 0.2 0.94

RCond = 0.182788206403613

Info = 0