Chapter 1: Linear Systems

lin_sol_gen_band

Solves a real general band system of linear equations, Ax = b. Using optional arguments,
any of sev­eral related computations can be performed. These extra tasks include computing
the LU factoriza­tion of A using partial pivoting, solving ATx = b, or computing the solution
of  Ax = b given the LU factorization of A.

Synopsis

#include <imsl.h>

float *imsl_f_lin_sol_gen_band (int n, float a[], int nlca, int nuca,
float b[], …, 0)

The type double procedure is imsl_d_lin_sol_gen_band.

Required Arguments

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

float a[]   (Input)
Array of size (nlca + nuca + 1) containing the n × n banded coefficient matrix in band storage mode.

int nlca   (Input)
Number of lower codiagonals in a.

int nuca   (Input)
Number of upper codiagonals in a.

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 free. If no solution was computed, then NULL is returned.

Synopsis with Optional Arguments

#include <imsl.h>

float *imsl_f_lin_sol_gen_band (int n, float a[], int nlca,
int nuca, float b[],
IMSL_TRANSPOSE,
IMSL_RETURN_USER, float x[],
IMSL_FACTOR, int **p_pvt, float **p_factor,
IMSL_FACTOR_USER, int pvt[], float factor[],
IMSL_CONDITION, float *condition,
IMSL_FACTOR_ONLY,
IMSL_SOLVE_ONLY,
IMSL_BLOCKING_FACTOR, int block_factor,
0)

Optional Arguments

IMSL_TRANSPOSE
Solve ATx = b.
Default: Solve Ax = b.

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

IMSL_FACTOR, int **p_pvt, float **p_factor   (Output)

p_pvt: The address of a pointer to an array of length n containing the pivot sequence for the factorization. On return, the necessary space is allocated by imsl_f_lin_sol_gen_band. Typically, int *p_pvt is declared and &p_pvt is used as an argument.

p_factor: The address of a pointer to an array of size
(2nlca + nuca + 1) × n containing the LU factorization of A with column pivoting. On return, the necessary space is allocated by imsl_f_lin_sol_gen_band. Typically, float *p_factor is declared and &p_factor is used as an argument.

IMSL_FACTOR_USER, int pvt[], float factor[]   (Input/Output)

pvt[]: A user-allocated array of size n containing the pivot sequence for the factorization.

factor[]: A user-allocated array of size (2nlca + nuca + 1) × n containing the LU factorization of A. The strictly lower triangular part of this array contains information necessary to construct L, and the upper triangular part contains U. If A is not needed, factor and a can share the first
(nlca + nuca + 1) × n locations.

These parameters are “Input” if IMSL_SOLVE_ONLY 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 LU factorization of A with partial pivoting. 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_gen_band is NULL.

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

IMSL_BLOCKING_FACTOR, int block_factor   (Input)
The blocking factor. block_factor must be set no larger than 32.
Default: block_factor = 1

Description

The function imsl_f_lin_sol_gen_band solves a system of linear algebraic equations with a real band matrix A. It first computes the LU factorization of A based on the blocked LU fac­torization algorithm given in Du Croz et al. (1990). Level-3 BLAS invocations are replaced with inline loops. The blocking factor block_factor has the default value of 1, but can be reset to any positive value not exceeding 32.

The solution of the linear system is then found by solving two simpler systems,
y = L-1b and x = U -1y. 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 ma­chine precision), a warning message is issued. This indicates that very small changes in A may pro­duce large changes in the solution x. The function imsl_f_lin_sol_gen_band fails if
U, the upper triangular part of the factorization, has a zero diagonal element.

Examples

Example 1

This example demonstrates the simplest use of this function by solving a system of four linear equations. This is the simplest usage of the function. The equations are as follows:

2x1 − x2 = 3

−3x1 + x2 − 2x3 = 1

−x3 + 2x4 = 11

2x3 + x4 = −2

#include <imsl.h>


void main ()

{

        int         n = 4;

        int         nuca = 1;

        int         nlca = 1;

        float      *x;


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


        float a[] = {0.0, -1.0, -2.0, 2.0,

                2.0, 1.0, -1.0, 1.0,

                -3.0, 0.0, 2.0, 0.0};

        float b[] = {3.0, 1.0, 11.0, -2.0};


        x = imsl_f_lin_sol_gen_band (n, a, nlca, nuca, 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

         2           1          -3           4

Example 2

In this example, the problem Ax = b is solved using the data from the first example. This time, the factorizations are returned and the problem ATx = b is solved without recomputing LU.

#include <imsl.h>


void main ()

{

        int         n = 4;

        int         nuca = 1;

        int         nlca = 1;

        int        *pivot;

        float       x[4];

        float      *factor;


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


        float a[] = {0.0, -1.0, -2.0, 2.0,

                     2.0, 1.0, -1.0, 1.0,

                    -3.0, 0.0, 2.0, 0.0};

        float b[] = {3.0, 1.0, 11.0, -2.0};


                        /*  Solve Ax = b and return LU */


        imsl_f_lin_sol_gen_band (n, a, nlca, nuca, b,

                IMSL_FACTOR, &pivot, &factor,

                IMSL_RETURN_USER, x,

                0);


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


                        /*  Use precomputed LU to solve trans(A)x = b */

                        /*  The original matrix A is not needed */


        imsl_f_lin_sol_gen_band (n, (float*) 0, nlca, nuca, b,

                IMSL_FACTOR_USER, pivot, factor,

                IMSL_SOLVE_ONLY,

                IMSL_TRANSPOSE,

                IMSL_RETURN_USER, x,

                0);


        imsl_f_write_matrix ("Solution of trans(A)x = b", 1, n, x, 0);

}

Output

              Solution of Ax = b

         1           2           3           4

         2           1          -3           4

 

           Solution of trans(A)x = b

         1           2           3           4

        -6          -5          -1          -0

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_SINGULAR_MATRIX                         The input matrix is singular.

Visual Numerics, Inc.
Visual Numerics - Developers of IMSL and PV-WAVE
http://www.vni.com/
PHONE: 713.784.3131
FAX:713.781.9260