lin_sol_posdef_band (complex)

Solves a complex Hermitian positive definite system of linear equations Ax = b in band symmetric storage mode. Using optional arguments, any of several related computations can be performed. These extra tasks include computing the RHR Cholesky factorization of A, computing the solution of Ax = b given the Cholesky factorization of A, or estimating the L1 condition number of A.

Synopsis

#include <imsl.h>

f_complex *imsl_c_lin_sol_posdef_band (int n, f_complex a[], int ncoda, f_complex b[], …, 0)

The type double function is imsl_z_lin_sol_posdef_band.

Required Arguments

int n (Input)
Number of rows and columns in the matrix.

f_complex a[] (Input)
Array of size (ncoda + 1) × n containing the n × n positive definite band coefficient matrix in band symmetric storage mode.

int ncoda (Input)
Number of upper codiagonals of the matrix.

f_complex b[] (Input)
Array of size n containing the right-hand side.

Return Value

A pointer to the solution x of the linear system Ax = b. To release this space use imsl_free. If no solution was computed, then NULL is returned.

Synopsis with Optional Arguments

#include <imsl.h>

f_complex *imsl_c_lin_sol_posdef_band (int n, f_complex a[], int ncoda, f_complex b[],

IMSL_RETURN_USER, f_complex x[],

IMSL_FACTOR, f_complex **p_factor,

IMSL_FACTOR_USER, f_complex factor[],

IMSL_CONDITION, float *condition,

IMSL_FACTOR_ONLY,

IMSL_SOLVE_ONLY,

Optional Arguments

IMSL_RETURN_USER, f_complex x[] (Output)
A user-allocated array of length n containing the solution x.

IMSL_FACTOR, f_complex **p_factor (Output)
The address of a pointer to an array of size (ncoda + 1) × n containing the RHR factorization of A. On return, the necessary space is allocated by imsl_c_lin_sol_posdef_band. Typically, f_complex *p_factor is declared and &p_factor is used as an argument.

IMSL_FACTOR_USER, f_complex factor[] (Input/Output)
A user-allocated array of size (ncoda + 1) × n containing the RHR factorization of A in band symmetric form. If A is not needed, factor and a can share the same storage. These parameters are “Input” if IMSL_SOLVE is specified. They are “Output” otherwise.

IMSL_CONDITION, float *condition (Output)
A pointer to a scalar containing an estimate of the L1 norm condition number of the matrix A. This option cannot be used with the option IMSL_SOLVE_ONLY.

IMSL_FACTOR_ONLY
Compute the RHR factorization of A. If IMSL_FACTOR_ONLY is used, either IMSL_FACTOR or IMSL_FACTOR_USER is required. The argument b is then ignored, and the returned value of imsl_c_lin_sol_posdef_band is NULL.

IMSL_SOLVE_ONLY
Solve Ax = b given the RHR factorization previously computed by imsl_c_lin_sol_posdef_band. By default, the solution to Ax = b is pointed to by imsl_c_lin_sol_posdef_band. If IMSL_SOLVE_ONLY is used, argument IMSL_FACTOR_USER is required and the argument a is ignored.

Description

The function imsl_c_lin_sol_posdef_band solves a system of linear algebraic equations with a real symmetric positive definite band coefficient matrix A. It computes the RHR Cholesky factorization of A. Argument R is an upper triangular band matrix.

When the solution to the linear system or the inverse of the matrix is sought, an estimate of the L1 condition number of A is computed using Higham’s modifications to Hager’s method, as given in Higham (1988). If the estimated condition number is greater than 1/ɛ (where ɛ is the machine precision), a warning message is issued. This indicates that very small changes in A may produce large changes in the solution x.

The function imsl_c_lin_sol_posdef_band fails if any submatrix of R is not positive definite or if R has a zero diagonal element. These errors occur only if A is very close to a singular matrix or to a matrix which is not positive definite.

The function imsl_c_lin_sol_posdef_band is based partially on the LINPACK subroutines SPBFA and CPBSL; see Dongarra et al. (1979).

Examples

Example 1

Solve a linear system Ax = b where

#include <imsl.h>

int main()

{

int n = 5;

int ncoda = 1;

f_complex *x;

/* Note that a is in band storage mode */

f_complex a[] =

{{0.0, 0.0}, {-1.0, 1.0}, {1.0, 2.0}, {0.0, 4.0},

{1.0, 1.0},

{2.0, 0.0}, {4.0, 0.0}, {10.0, 0.0}, {6.0, 0.0},

{9.0, 0.0}};

f_complex b[] =

{{1.0, 5.0}, {12.0, -6.0}, {1.0, -16.0},{-3.0, -3.0},

{25.0, 16.0}};

x = imsl_c_lin_sol_posdef_band (n, a, ncoda, b, 0);

imsl_c_write_matrix ("Solution, x, of Ax = b", n, 1, x, 0);

}

Output

Solution, x, of Ax = b

1 ( 2, 1)

2 ( 3, -0)

3 ( -1, -1)

4 ( 0, -2)

5 ( 3, 2)

Example 2

This example solves the same problem Ax = b given in the first example. The solution is returned in user-allocated space and an estimate of κ1(A) is computed. Additionally, the RHR factorization is returned. Then, knowing that κ1(A) = ∥A∥ ∥A-1∥, the condition number is computed directly and compared to the estimate from Higham’s method.

#include <imsl.h>

#include <stdio.h>

#include <math.h>

int main()

{

int n = 5, ncoda = 1, i, j;

/* Note that a is in band storage mode */

f_complex a[] =

{{0.0, 0.0}, {-1.0, 1.0}, {1.0, 2.0}, {0.0, 4.0},

{1.0, 1.0},

{2.0, 0.0}, {4.0, 0.0}, {10.0, 0.0}, {6.0, 0.0},

{9.0, 0.0}};

f_complex b[] =

{{1.0, 5.0}, {12.0, -6.0}, {1.0, -16.0},{-3.0, -3.0},

{25.0, 16.0}};

f_complex x[5], e_i[5], *factor;

float condition, column_norm, inverse_norm;

imsl_c_lin_sol_posdef_band (n, a, ncoda, b,

IMSL_FACTOR, &factor,

IMSL_CONDITION, &condition,

IMSL_RETURN_USER, x,

0);

imsl_c_write_matrix ("Solution, x, of Ax = b", 1, n, x, 0);

/* Find one norm of inverse */

inverse_norm = 0.0;

for (i=0; i<n; i++) {

for (j=0; j<n; j++) e_i[j] = imsl_cf_convert (0.0, 0.0);

e_i[i] = imsl_cf_convert (1.0, 0.0);

/* Determine one norm of each column of inverse */

imsl_c_lin_sol_posdef_band (n, a, ncoda, e_i,

IMSL_FACTOR_USER, factor,

IMSL_SOLVE_ONLY,

IMSL_RETURN_USER, x,

0);

column_norm = imsl_c_vector_norm (n, x,

IMSL_ONE_NORM,

0);

/* The max of the column norms is the norm of inv(A) */

if (inverse_norm < column_norm)

inverse_norm = column_norm;

}

/* By observation, one norm of A is 14+sqrt(5) */

printf ("\nHigham’s condition estimate = %7.4f\n", condition);

printf ("Direct condition estimate = %7.4f\n",

(14.0+sqrt(5.0))*inverse_norm);

}

Output

Solution, x, of Ax = b

1 2 3

( 2, 1) ( 3, -0) ( -1, -1)

4 5

( 0, -2) ( 3, 2)

Higham’s condition estimate = 19.3777

Direct condition estimate = 19.3777

Warning Errors

IMSL_ILL_CONDITIONED

The input matrix is too ill-conditioned. An estimate of the reciprocal of its L1 condition number is “rcond” = #.
The solution might not be accurate.

Fatal Errors

IMSL_NONPOSITIVE_MATRIX	The leading # by # submatrix of the input matrix is not positive definite.
IMSL_SINGULAR_MATRIX	The input matrix is singular.