Chapter 1: Linear Systems

lin_sol_posdef

Solves a real symmetric positive definite system of linear equations Ax = b. Using optional arguments, any of several related computations can be performed. These extra tasks include computing the Cholesky factor, L, of A such that A = LLT, computing the inverse matrix A-1, or computing the solution of Ax = b given the Cholesky factor, L.

Synopsis

#include <imsl.h>

float *imsl_f_lin_sol_posdef (int n, float a[], float b[], …, 0)

The type double procedure is imsl_d_lin_sol_posdef.

Required Arguments

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

float a[]   (Input)
Array of size n × n containing the matrix.

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

Return Value

A pointer to the solution x of the symmetric positive definite 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_posdef (int n, float a[], float b[],
IMSL_A_COL_DIM, int a_col_dim,
IMSL_RETURN_USER, float x[],
IMSL_FACTOR, float **p_factor,
IMSL_FACTOR_USER, float factor[],
IMSL_FAC_COL_DIM, int fac_col_dim,
IMSL_INVERSE, float **p_inva,
IMSL_INVERSE_USER, float inva[],
IMSL_INV_COL_DIM, int inv_col_dim,
IMSL_CONDITION, float *cond,
IMSL_FACTOR_ONLY,
IMSL_SOLVE_ONLY,
IMSL_INVERSE_ONLY,
0)

Optional Arguments

IMSL_A_COL_DIM, int a_col_dim   (Input)
The column dimension of the array a.
Default: a_col_dim = n

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 n × n containing the
LLT factorization of A. On return, the necessary space is allocated by imsl_f_lin_sol_posdef. The lower-triangular part of this array contains
L and the upper-triangular part contains LT. 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 n × n containing the LLT factorization of A.
The lower-triangular part of this array contains L, and the upper-triangular part contains LT. If A is not needed, a and factor can share the same storage.
If IMSL_SOLVE is specified, it is input; otherwise, it is output.

IMSL_FAC_COL_DIM, int fac_col_dim   (Input)
The column dimension of the array containing the LLT factorization of A.
Default: fac_col_dim = n

IMSL_INVERSE, float **p_inva   (Output)
The address of a pointer to an array of size n × n containing the inverse of the matrix A. On return, the necessary space is allocated by imsl_f_lin_sol_posdef. Typically,  float *p_inva is declared, and &p_inva is used as an argument.

IMSL_INVERSE_USER, float inva[]   (Output)
A user-allocated array of size n × n containing the inverse of A.

IMSL_INV_COL_DIM, int inva_col_dim   (Input)
The column dimension of the array containing the inverse of A.
Default: inva_col_dim = n

IMSL_CONDITION, float *cond   (Output)
A pointer to a scalar containing an estimate of the L1 norm condition number of the matrix A. Do not use this option with IMSL_SOLVE_ONLY.

IMSL_FACTOR_ONLY
Compute the Cholesky factorization LLT 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 is NULL.

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

IMSL_INVERSE_ONLY
Compute the inverse of the matrix A. If IMSL_INVERSE_ONLY is used, either IMSL_INVERSE or IMSL_INVERSE_USER is required. The argument b is then ignored, and the returned value of imsl_f_lin_sol_posdef is NULL.

Description

The function imsl_f_lin_sol_posdef solves a system of linear algebraic equations having a symmetric positive definite coefficient matrix A. The function first computes the Cholesky factorization LLT of A. The solution of the linear system is then found by solving the two simpler systems, y = L-1b and x = L-Ty. 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 the same algorithm as in Dongarra et al. (1979). 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 fails if L, the lower-triangular matrix in the factorization, has a zero diagonal element.

Examples

Example 1

A system of three linear equations with a symmetric positive definite coefficient matrix is solved in this example. The equations are listed below:

x1 − 3x2 + 2x3 = 27

−3x1 + 10x2 − 5x3 = −78

2x1 − 5x2 + 6x3 = 64

#include <imsl.h>

main()
{
    int         n = 3;
    float       *x;
    float       a[] =  {1.0, -3.0,  2.0,
                        -3.0, 10.0, -5.0,
                         2.0, -5.0,  6.0};
    float       b[] =  {27.0, -78.0, 64.0};

                                /* Solve Ax = b  for  x  */
    x = imsl_f_lin_sol_posdef (n, a, b, 0);

                                /* Print x */
    imsl_f_write_matrix ("Solution, x, of Ax = b", 1, n, x, 0);
}

Output

Solution, x, of Ax = b
   1           2           3
   1          -4           7

Example 2

This example solves the same system of three linear equations as in the initial example, but this time returns the LLT factorization of A. The solution x is returned in an array allocated in the main program.

#include <imsl.h>

main()
{
    int         n = 3;
    float       x[3], *p_factor;
    float       a[] =  {1.0, -3.0,  2.0,
                        -3.0, 10.0, -5.0,
                         2.0, -5.0,  6.0};
    float       b[] =  {27.0, -78.0, 64.0};

                                /* Solve Ax = b  for  x  */
    imsl_f_lin_sol_posdef (n, a, b,
                IMSL_RETURN_USER, x,
                IMSL_FACTOR, &p_factor,
                0);

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

                                /* Print Cholesky factor of A */
    imsl_f_write_matrix ("Cholesky factor L, and trans(L), of A",
                         n, n, p_factor, 0);
}

Output

Solution, x, of Ax = b
1           2           3
1          -4           7
 
Cholesky factor L, and trans(L), of A
            1           2           3
1           1          -3           2
2          -3           1           1
3           2           1           1

Example 3

This example solves the same system as in the initial example, but given the Cholesky factors of A.

#include <imsl.h>

main()
{
    int         n = 3;
    float       *x, *a;
    float       factor[ ] =  {1.0, -3.0, 2.0,
                              -3.0, 1.0, 1.0,
                              2.0,  1.0, 1.0};
    float       b[ ] =  {27.0, -78.0, 64.0};

                                /* Solve Ax = b  for  x  */
    x = imsl_f_lin_sol_posdef (n, a, b,
                    IMSL_FACTOR_USER, factor,
                    IMSL_SOLVE_ONLY,
                    0);

                                /* Print x */
    imsl_f_write_matrix ("Solution, x, of Ax = b", 1, n, x, 0);
}

Output

Solution, x, of Ax = b
1           2           3
1          -4           7

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.

IMSL_SINGULAR_TRI_MATRIX                The input triangular matrix is singular. The index of the first zero diagonal element is #.

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