lin_sol_posdef_band
Solves a real symmetric 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 RTR 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>
float *imsl_f_lin_sol_posdef_band (int n, float a[], int ncoda, float b[], …, 0)
The type double function is imsl_d_lin_sol_posdef_band.
Required Arguments
int n (Input)
Number of rows and columns in the matrix.
float 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.
float 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>
float *imsl_f_lin_sol_posdef_band (int n, float a[], int ncoda, float b[],
IMSL_RETURN_USER, float x[],
IMSL_FACTOR, float **p_factor,
IMSL_FACTOR_USER, float factor[],
IMSL_CONDITION, float *cond,
IMSL_FACTOR_ONLY,
IMSL_SOLVE_ONLY,
0)
Optional Arguments
IMSL_RETURN_USER, float x[] (Output)
A user-allocated array of length n containing the solution x.
IMSL_FACTOR, float **p_factor (Output)
The address of a pointer to an array of size (ncoda + 1) × n containing the LLT factorization of A. On return, the necessary space is allocated by imsl_f_lin_sol_posdef_band. Typically, float *p_factor is declared and &p_factor is used as an argument.
IMSL_FACTOR_USER, float factor[] (Input/Output)
A user-allocated array of size (ncoda + 1) × n containing the LLT 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 *cond (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 LLT 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_f_lin_sol_posdef_band is NULL.
IMSL_SOLVE_ONLY
Solve Ax = b given the LLT factorization previously computed by imsl_f_lin_sol_posdef_band. By default, the solution to Ax = b is pointed to by imsl_f_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_f_lin_sol_posdef_band solves a system of linear algebraic equations with a real symmetric positive definite band coefficient matrix A. It computes the RTR Cholesky factorization of A. 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_f_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_f_lin_sol_posdef_band is partially based on the LINPACK subroutines CPBFA and SPBSL; see Dongarra et al. (1979).
Example 1
Solves a system of linear equations Ax = b, where
#include <imsl.h>
int main()
{
int n = 4;
int ncoda = 2;
float *x;
/* Note that a is in band storage mode */
float a[] = {0.0, 0.0, -1.0, 1.0,
0.0, 0.0, 2.0, -1.0,
2.0, 4.0, 7.0, 3.0};
float b[] = {6.0, -11.0, -11.0, 19.0};
x = imsl_f_lin_sol_posdef_band (n, a, ncoda, b, 0);
imsl_f_write_matrix ("Solution, x, of Ax = b", 1, n, x, 0);
}
Output
Solution, x, of Ax = b
1 2 3 4
4 -6 2 9
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 RTR 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>
int main()
{
int n = 4;
int ncoda = 2;
float a[] =
{0.0, 0.0, -1.0, 1.0,
0.0, 0.0, 2.0, -1.0,
2.0, 4.0, 7.0, 3.0};
float b[] = {6.0, -11.0, -11.0, 19.0};
float x[4];
float e_i[4];
float *factor;
float condition;
float column_norm;
float inverse_norm;
int i;
int j;
imsl_f_lin_sol_posdef_band (n, a, ncoda, b,
IMSL_FACTOR, &factor,
IMSL_CONDITION, &condition,
IMSL_RETURN_USER, x,
0);
imsl_f_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] = 0.0;
e_i[i] = 1.0;
/* determine one norm of each column of inverse */
imsl_f_lin_sol_posdef_band (n, a, ncoda, e_i,
IMSL_FACTOR_USER, factor,
IMSL_SOLVE_ONLY,
IMSL_RETURN_USER, x,
0);
column_norm = imsl_f_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 11 */
printf ("\nHigham’s condition estimate = %f\n", condition);
printf ("Direct condition estimate = %f\n",
11.0*inverse_norm);
}
Output
Solution, x, of Ax = b
1 2 3 4
4 -6 2 9
Higham’s condition estimate = 8.650485
Direct condition estimate = 8.650485
Warning Errors
|
IMSL_ILL_CONDITIONED |
The input matrix is too ill-conditioned. An estimate of the reciprocal of its L1 condition number is “rcond” = #. |
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. |
