Integer type: int32 int64 nag_int show int32 show int32 show int64 show int64 show nag_int show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)

Chapter Contents

Chapter Introduction

NAG Toolbox

NAG Toolbox: nag_lapack_zpbsvx (f07hp)

▸▿ Contents

1 Purpose

2 Syntax

3 Description

4 References

▸▿ 5 Parameters

5.1 Compulsory Input Parameters

5.2 Optional Input Parameters

5.3 Output Parameters

6 Error Indicators and Warnings

7 Accuracy

8 Further Comments

9 Example

Purpose

nag_lapack_zpbsvx (f07hp) uses the Cholesky factorization

A = U^{H} U or A = L L^{H}

to compute the solution to a complex system of linear equations

A X = B,

where

A

is an

n

n

Hermitian positive definite band matrix of bandwidth

(2 k_{d} + 1)

and

X

and

B

are

n

r

matrices. Error bounds on the solution and a condition estimate are also provided.

Syntax

[ab, afb, equed, s, b, x, rcond, ferr, berr, info] = f07hp(fact, uplo, kd, ab, afb, equed, s, b, 'n', n, 'nrhs_p', nrhs_p)

[ab, afb, equed, s, b, x, rcond, ferr, berr, info] = nag_lapack_zpbsvx(fact, uplo, kd, ab, afb, equed, s, b, 'n', n, 'nrhs_p', nrhs_p)

Description

nag_lapack_zpbsvx (f07hp) performs the following steps:

fact ='E'

, real diagonal scaling factors,

D_{S}

, are computed to equilibrate the system:

(D_{S} A D_{S}) (D_{S}^{- 1} X) = D_{S} B .

Whether or not the system will be equilibrated depends on the scaling of the matrix

A

, but if equilibration is used,

A

is overwritten by

D_{S} A D_{S}

and

B

D_{S} B

fact ='N'

'E'

, the Cholesky decomposition is used to factor the matrix

A

(after equilibration if

fact ='E'

) as

A = U^{H} U

uplo ='U'

A = L L^{H}

uplo ='L'

, where

U

is an upper triangular matrix and

L

is a lower triangular matrix.

If the leading

i

i

principal minor of

A

is not positive definite, then the function returns with

info = i

. Otherwise, the factored form of

A

is used to estimate the condition number of the matrix

A

. If the reciprocal of the condition number is less than machine precision,

info \geq n + 1

is returned as a warning, but the function still goes on to solve for

X

and compute error bounds as described below.

The system of equations is solved for

X

using the factored form of

A

Iterative refinement is applied to improve the computed solution matrix and to calculate error bounds and backward error estimates for it.

If equilibration was used, the matrix

X

is premultiplied by

D_{S}

so that it solves the original system before equilibration.

References

Anderson E, Bai Z, Bischof C, Blackford S, Demmel J, Dongarra J J, Du Croz J J, Greenbaum A, Hammarling S, McKenney A and Sorensen D (1999) LAPACK Users' Guide (3rd Edition) SIAM, Philadelphia http://www.netlib.org/lapack/lug

Golub G H and Van Loan C F (1996) Matrix Computations (3rd Edition) Johns Hopkins University Press, Baltimore

Higham N J (2002) Accuracy and Stability of Numerical Algorithms (2nd Edition) SIAM, Philadelphia

Parameters

Compulsory Input Parameters

1: $fact$ – string (length ≥ 1)

Specifies whether or not the factorized form of the matrix

A

is supplied on entry, and if not, whether the matrix

A

should be equilibrated before it is factorized.

$fact ='F'$: afb contains the factorized form of $A$ . If $equed ='Y'$ , the matrix $A$ has been equilibrated with scaling factors given by s. ab and afb will not be modified.
$fact ='N'$: The matrix $A$ will be copied to afb and factorized.
$fact ='E'$: The matrix $A$ will be equilibrated if necessary, then copied to afb and factorized.

Constraint:

fact ='F'

'N'

'E'

2: $uplo$ – string (length ≥ 1)

uplo ='U'

, the upper triangle of

A

is stored.

uplo ='L'

, the lower triangle of

A

is stored.

Constraint:

uplo ='U'

'L'

3: $kd$ – int64int32nag_int scalar

k_{d}

, the number of superdiagonals of the matrix

A

uplo ='U'

, or the number of subdiagonals if

uplo ='L'

Constraint:

kd \geq 0

4: $ab (ldab :)$ – complex array

The first dimension of the array ab must be at least

kd + 1

The second dimension of the array ab must be at least

\max (1, n)

The upper or lower triangle of the Hermitian band matrix

A

, except if

fact ='F'

and

equed ='Y'

, in which case ab must contain the equilibrated matrix

D_{S} A D_{S}

The matrix is stored in rows

1

k_{d} + 1

, more precisely,

if $uplo ='U'$ , the elements of the upper triangle of $A$ within the band must be stored with element $A_{i j}$ in $ab (k_{d} + 1 + i - j, j) for \max (1, j - k_{d}) \leq i \leq j$ ;
if $uplo ='L'$ , the elements of the lower triangle of $A$ within the band must be stored with element $A_{i j}$ in $ab (1 + i - j, j) for j \leq i \leq \min (n, j + k_{d}) .$

5: $afb (ldafb :)$ – complex array

The first dimension of the array afb must be at least

kd + 1

The second dimension of the array afb must be at least

\max (1, n)

fact ='F'

, afb contains the triangular factor

U

L

from the Cholesky factorization

A = U^{H} U

A = L L^{H}

of the band matrix

A

, in the same storage format as

A

. If

equed ='Y'

, afb is the factorized form of the equilibrated matrix

A

6: $equed$ – string (length ≥ 1)

fact ='N'

'E'

, equed need not be set.

fact ='F'

, equed must specify the form of the equilibration that was performed as follows:

if $equed ='N'$ , no equilibration;
if $equed ='Y'$ , equilibration was performed, i.e., $A$ has been replaced by $D_{S} A D_{S}$ .

Constraint: if

fact ='F'

equed ='N'

'Y'

7: $s (:)$ – double array

The dimension of the array s must be at least

\max (1, n)

fact ='N'

'E'

, s need not be set.

fact ='F'

and

equed ='Y'

, s must contain the scale factors,

D_{S}

, for

A

; each element of s must be positive.

8: $b (ldb :)$ – complex array

The first dimension of the array b must be at least

\max (1, n)

The second dimension of the array b must be at least

\max (1, nrhs_p)

The

n

r

right-hand side matrix

B

Optional Input Parameters

1: $n$ – int64int32nag_int scalar: Default: the first dimension of the array b and the second dimension of the arrays ab, afb, s.
$n$ , the number of linear equations, i.e., the order of the matrix $A$ .

Constraint: $n \geq 0$ .
2: $nrhs_p$ – int64int32nag_int scalar: Default: the second dimension of the array b.
$r$ , the number of right-hand sides, i.e., the number of columns of the matrix $B$ .

Constraint: $nrhs_p \geq 0$ .

Output Parameters

1: $ab (ldab :)$ – complex array: The first dimension of the array ab will be $kd + 1$ .
The second dimension of the array ab will be $\max (1, n)$ .

If $fact ='E'$ and $equed ='Y'$ , ab stores $D_{S} A D_{S}$ .
2: $afb (ldafb :)$ – complex array: The first dimension of the array afb will be $kd + 1$ .
The second dimension of the array afb will be $\max (1, n)$ .

If $fact ='N'$ , afb returns the triangular factor $U$ or $L$ from the Cholesky factorization $A = U^{H} U$ or $A = L L^{H}$ .
If $fact ='E'$ , afb returns the triangular factor $U$ or $L$ from the Cholesky factorization $A = U^{H} U$ or $A = L L^{H}$ of the equilibrated matrix $A$ (see the description of ab for the form of the equilibrated matrix).
3: $equed$ – string (length ≥ 1): If $fact ='F'$ , equed is unchanged from entry.
Otherwise, if no constraints are violated, equed specifies the form of the equilibration that was performed as specified above.
4: $s (:)$ – double array: The dimension of the array s will be $\max (1, n)$

If $fact ='F'$ , s is unchanged from entry.
Otherwise, if no constraints are violated and $equed ='Y'$ , s contains the scale factors, $D_{S}$ , for $A$ ; each element of s is positive.
5: $b (ldb :)$ – complex array: The first dimension of the array b will be $\max (1, n)$ .
The second dimension of the array b will be $\max (1, nrhs_p)$ .

If $equed ='N'$ , b is not modified.
If $equed ='Y'$ , b stores $D_{S} B$ .
6: $x (ldx :)$ – complex array: The first dimension of the array x will be $\max (1, n)$ .
The second dimension of the array x will be $\max (1, nrhs_p)$ .

If $info = 0$ or $n + 1$ , the $n$ by $r$ solution matrix $X$ to the original system of equations. Note that the arrays $A$ and $B$ are modified on exit if $equed ='Y'$ , and the solution to the equilibrated system is $D_{S}^{- 1} X$ .
7: $rcond$ – double scalar: If no constraints are violated, an estimate of the reciprocal condition number of the matrix $A$ (after equilibration if that is performed), computed as $rcond = 1.0 / ({‖A‖}_{1} {‖A^{- 1}‖}_{1})$ .
8: $ferr (nrhs_p)$ – double array: If $info = 0$ or $n + 1$ , an estimate of the forward error bound for each computed solution vector, such that ${‖{\hat{x}}_{j} - x_{j}‖}_{\infty} / {‖x_{j}‖}_{\infty} \leq ferr (j)$ where ${\hat{x}}_{j}$ is the $j$ th column of the computed solution returned in the array x and $x_{j}$ is the corresponding column of the exact solution $X$ . The estimate is as reliable as the estimate for rcond, and is almost always a slight overestimate of the true error.
9: $berr (nrhs_p)$ – double array: If $info = 0$ or $n + 1$ , an estimate of the component-wise relative backward error of each computed solution vector ${\hat{x}}_{j}$ (i.e., the smallest relative change in any element of $A$ or $B$ that makes ${\hat{x}}_{j}$ an exact solution).
10: $info$ – int64int32nag_int scalar: $info = 0$ unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

$info < 0$: If $info = - i$ , argument $i$ had an illegal value. An explanatory message is output, and execution of the program is terminated.

$info > 0 and info \leq n$: The leading minor of order $_$ of $A$ is not positive definite, so the factorization could not be completed, and the solution has not been computed. $rcond = 0.0$ is returned.

W $info = n + 1$: $U$ (or $L$ ) is nonsingular, but rcond is less than machine precision, meaning that the matrix is singular to working precision. Nevertheless, the solution and error bounds are computed because there are a number of situations where the computed solution can be more accurate than the value of rcond would suggest.

Accuracy

For each right-hand side vector

b

, the computed solution

x

is the exact solution of a perturbed system of equations

(A + E) x = b

, where

if $uplo ='U'$ , $|E| \leq c (n) ε |U^{H}| |U|$ ;
if $uplo ='L'$ , $|E| \leq c (n) ε |L| |L^{H}|$ ,

c (n)

is a modest linear function of

n

, and

ε

is the machine precision. See Section 10.1 of Higham (2002) for further details.

\hat{x}

is the true solution, then the computed solution

x

satisfies a forward error bound of the form

\frac{{‖x - \hat{x}‖}_{\infty}}{{‖\hat{x}‖}_{\infty}} \leq w_{c} cond (A, \hat{x}, b)

where

cond (A, \hat{x}, b) = {‖|A^{- 1}| (|A| |\hat{x}| + |b|)‖}_{\infty} / {‖\hat{x}‖}_{\infty} \leq cond (A) = {‖|A^{- 1}| |A|‖}_{\infty} \leq κ_{\infty} (A)

. If

\hat{x}

is the

j

th column of

X

, then

w_{c}

is returned in

berr (j)

and a bound on

{‖x - \hat{x}‖}_{\infty} / {‖\hat{x}‖}_{\infty}

is returned in

ferr (j)

. See Section 4.4 of Anderson et al. (1999) for further details.

Further Comments

When

n ≫ k

, the factorization of

A

requires approximately

4 n {(k + 1)}^{2}

floating-point operations, where

k

is the number of superdiagonals.

For each right-hand side, computation of the backward error involves a minimum of

32 n k

floating-point operations. Each step of iterative refinement involves an additional

48 n k

operations. At most five steps of iterative refinement are performed, but usually only one or two steps are required. Estimating the forward error involves solving a number of systems of equations of the form

A x = b

; the number is usually

4

5

and never more than

11

. Each solution involves approximately

16 n k

operations.

The real analogue of this function is nag_lapack_dpbsvx (f07hb).

Example

This example solves the equations

A X = B,

where

A

is the Hermitian positive definite band matrix

A = (\begin{array}{r} 9.39 & 1.08 - 1.73 i & 0 & 0 \\ 1.08 + 1.73 i & 1.69 & - 0.04 + 0.29 i & 0 \\ 0 & - 0.04 - 0.29 i & 2.65 & - 0.33 + 2.24 i \\ 0 & 0 & - 0.33 - 2.24 i & 2.17 \end{array})

and

B = (\begin{array}{r} - 12.42 + 68.42 i & 54.30 - 56.56 i \\ - 9.93 + 0.88 i & 18.32 + 4.76 i \\ - 27.30 - 0.01 i & - 4.40 + 9.97 i \\ 5.31 + 23.63 i & 9.43 + 1.41 i \end{array}) .

Error estimates for the solutions, information on equilibration and an estimate of the reciprocal of the condition number of the scaled matrix

A

are also output.

Open in the MATLAB editor: f07hp_example

function f07hp_example


fprintf('f07hp example results\n\n');

uplo = 'U';
kd = int64(1);
n  = int64(4);
ab = [0,          1.08 - 1.73i,  -0.04 + 0.29i,  -0.33 + 2.24i;
      9.39 + 0i,  1.69 + 0i,      2.65 + 0i,      2.17 + 0i];

% Input parameters
fact  = 'Equilibration';
afb   = complex(zeros(kd+1, n));
equed = ' ';
s     = zeros(n,1);

% RHS
b = [-12.42 + 68.42i,  54.30 - 56.56i;
      -9.93 +  0.88i,  18.32 +  4.76i;
     -27.30 -  0.01i,  -4.40 +  9.97i;
       5.31 + 23.63i,   9.43 +  1.41i];

% Solve
[ab, afb, equed, s, b, x, rcond, ferr, berr, info] = ...
  f07hp( ...
         fact, uplo, kd, ab, afb, equed, s, b);

disp('Solution(s)');
disp(x);
disp('Backward errors (machine-dependent)');
fprintf('%10.1e',berr);
fprintf('\n');
disp('Estimated forward error bounds (machine-dependent)');
fprintf('%10.1e',ferr);
fprintf('\n\n');
disp('Estimate of reciprocal condition number');
fprintf('%10.1e\n\n',rcond);
if equed=='N'
  fprintf('A has not been equilibrated\n');
else
  fprintf('A has been equilibrated\n');
end

f07hp example results

Solution(s)
  -1.0000 + 8.0000i   5.0000 - 6.0000i
   2.0000 - 3.0000i   2.0000 + 3.0000i
  -4.0000 - 5.0000i  -8.0000 + 4.0000i
   7.0000 + 6.0000i  -1.0000 - 7.0000i

Backward errors (machine-dependent)
   9.1e-17   9.0e-17
Estimated forward error bounds (machine-dependent)
   3.7e-14   3.2e-14

Estimate of reciprocal condition number
   7.6e-03

A has not been equilibrated

PDF version (NAG web site, 64-bit version, 64-bit version)

Chapter Contents

Chapter Introduction

NAG Toolbox