lin_sol_nonnegdef
Solves a real symmetric nonnegative definite system of linear equations Ax = b. Using options, computes a Cholesky factorization of the matrix A, such that A = RTR = LLT. Computes the solution to Ax = b given the Cholesky factor.
Synopsis
#include <imsl.h>
float *imsl_f_lin_sol_nonnegdef (int n, float a[], float b[], …, 0)
The type double function is imsl_d_lin_sol_nonnegdef.
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
Using required arguments, imsl_f_lin_sol_nonnegdef returns a pointer to a solution x of the linear system. To release this space, use imsl_free. If no value can be computed, NULL is returned.
Synopsis with Optional Arguments
#include <imsl.h>
float *imsl_f_lin_sol_nonnegdef (int n, float a[], float b[],
IMSL_RETURN_USER, float x[],
IMSL_A_COL_DIM, int a_col_dim,
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_TOLERANCE, float tol,
IMSL_FACTOR_ONLY,
IMSL_SOLVE_ONLY,
IMSL_INVERSE_ONLY,
0)
Optional Arguments
IMSL_RETURN_USER, float x[] (Output)
A user-allocated array of length n containing the solution x. When this option is specified, no storage is allocated for the solution, and imsl_f_lin_sol_nonnegdef returns a pointer to the array x.
IMSL_A_COL_DIM, int a_col_dim (Input)
The column dimension of the array a.
Default: a_col_dim = n
IMSL_FACTOR, float **p_factor (Output)
The address of a pointer to an array of size n × n containing the LLT factorization of A. When this option is specified, the space for the factor matrix is allocated by imsl_f_lin_sol_nonnegdef. The lower-triangular part of the factor array contains L, and the upper-triangular part contains LTR. 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 factor contains L, and the upper-triangular part contains LT. If a is not needed, a and factor can be the same storage locations. If IMSL_SOLVE is specified, this parameter 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.
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 A. The space for this array is allocated by imsl_f_lin_sol_nonnegdef. 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. If a is not needed, a and factor can be the same storage locations. The storage locations for A cannot be the factorization and the inverse of A at the same time.
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_TOLERANCE,
float tol (Input)
Tolerance used in determining linear dependence. See the documentation for
-imsl_f_machine (
imsl_f_machine(float)) in Chapter 12, “Utilities.”
Default:
tol = 100
* imsl_f_machine(4) IMSL_FACTOR_ONLY
Compute the LLT factorization of A only. The argument b is ignored, and either the optional argument IMSL_FACTOR or IMSL_FACTOR_USER is required.
IMSL_SOLVE_ONLY
Solve Ax = b using the factorization previously computed by this function. The argument a is ignored, and the optional argument IMSL_FACTOR_USER is required.
IMSL_INVERSE_ONLY
Compute the inverse of A only. The argument b is ignored, and either the optional argument IMSL_INVERSE or IMSL_INVERSE_USER is required.
Description
The function imsl_f_lin_sol_nonnegdef solves a system of linear algebraic equations having a symmetric nonnegative definite (positive semidefinite) coefficient matrix. It first computes a Cholesky (LLT or RTR) factorization of the coefficient matrix A.
The factorization algorithm is based on the work of Healy (1968) and proceeds sequentially by columns. The i-th column is declared to be linearly dependent on the first i − 1 columns if
where ε (specified in tol) may be set by the user. When a linear dependence is declared, all elements in the i-th row of R (column of L) are set to zero.
Modifications due to Farebrother and Berry (1974) and Barrett and Healy (1978) for checking for matrices that are not nonnegative definite also are incorporated. The function imsl_f_lin_sol_nonnegdef declares A to not be nonnegative definite and issues an error message if either of the following conditions are satisfied:
Healy’s (1968) algorithm and the function imsl_f_lin_sol_nonnegdef permit the matrices A and R to occupy the same storage. Barrett and Healy (1978) in their remark neglect this fact. The function imsl_f_lin_sol_nonnegdef uses
for aii in the above condition 2 to remedy this problem.
If an inverse of the matrix A is required and the matrix is not (numerically) positive definite, then the resulting inverse is a symmetric g2 inverse of A. For a matrix G to be a g2 inverse of a matrix A, G must satisfy conditions 1 and 2 for the Moore-Penrose inverse, but generally fail conditions 3 and 4. The four conditions for G to be a Moore-Penrose inverse of A are as follows:
1. AGA = A
2. GAG = G
3. AG is symmetric
4. GA is symmetric
The solution of the linear system Ax = b is computed by solving the factored version of the linear system RTRx = b as two successive triangular linear systems. In solving the triangular linear systems, if the elements of a row of R are all zero, the corresponding element of the solution vector is set to zero. For a detailed description of the algorithm, see Section 2 in Sallas and Lionti (1988).
Examples
Example 1
A solution to a system of four linear equations is obtained. Maindonald (1984, pp. 83-86 and 104-105) discusses the computations for the factorization and solution to this problem.
#include <imsl.h>
int main()
{
int n = 4;
float *x;
float a[] = {36.0, 12.0, 30.0, 6.0,
12.0, 20.0, 2.0, 10.0,
30.0, 2.0, 29.0, 1.0,
6.0, 10.0, 1.0, 14.0};
float b[] = {18.0, 22.0, 7.0, 20.0};
/* Solve Ax = b for x */
x = imsl_f_lin_sol_nonnegdef(n, a, b, 0);
/* Print solution, x, of Ax = b */
imsl_f_write_matrix("Solution, x", 1, n, x, 0);
}
Output
Solution, x
1 2 3 4
0.167 0.500 0.000 1.000
Example 2
The symmetric nonnegative definite matrix in the initial example is used to compute the factorization only in the first call to lin_sol_nonnegdef. The space needed for the factor is provided by the user. On the second call, both the LLT factorization and the right-hand side vector in the first example are used as the input to compute a solution x. It also illustrates another way to obtain the solution array x.
#include <imsl.h>
int main()
{
int n = 4, a_col_dim = 6;
float factor[36], x[5];
float a[] = {36.0, 12.0, 30.0, 6.0,
12.0, 20.0, 2.0, 10.0,
30.0, 2.0, 29.0, 1.0,
6.0, 10.0, 1.0, 14.0};
float b[] = {18.0, 22.0, 7.0, 20.0};
/* Factor A */
imsl_f_lin_sol_nonnegdef(n, a, b,
IMSL_FACTOR_USER, factor,
IMSL_FAC_COL_DIM, a_col_dim,
IMSL_FACTOR_ONLY,
0);
/* NULL is returned in */
/* this case. Another */
/* way to obtain the */
/* factor is to use the */
/* IMSL_FACTOR option. */
imsl_f_write_matrix("factor", n, n, factor,
IMSL_A_COL_DIM, a_col_dim,
0);
/* Get the solution using */
/* the factorized matrix. */
imsl_f_lin_sol_nonnegdef(n, a, b,
IMSL_FACTOR_USER, factor,
IMSL_FAC_COL_DIM, a_col_dim,
IMSL_RETURN_USER, x,
IMSL_SOLVE_ONLY,
0);
imsl_f_write_matrix("Solution, x, of Ax = b", 1, n, x, 0);
}
Output
Factor
1 2 3 4
1 6 2 5 1
2 2 4 -2 2
3 5 -2 0 0
4 1 2 0 3
Solution, x, of Ax = b
1 2 3 4
0.167 0.500 0.000 1.000
Example 3
This example uses the IMSL_INVERSE option to compute the symmetric g inverse of the symmetric nonnegative matrix in the first example. Maindonald (1984, p. 106) discusses the computations for this problem.
#include <imsl.h>
int main()
{
int n = 4;
float *p_a_inva, *p_a_inva_a, *p_inva;
float a[] =
{36.0, 12.0, 30.0, 6.0,
12.0, 20.0, 2.0, 10.0,
30.0, 2.0, 29.0, 1.0,
6.0, 10.0, 1.0, 14.0};
/* Get g2_inverse(a) */
imsl_f_lin_sol_nonnegdef(n, a, NULL,
IMSL_INVERSE, &p_inva,
IMSL_INVERSE_ONLY,
0);
/* Form a*g2_inverse(a) */
p_a_inva = imsl_f_mat_mul_rect("A*B",
IMSL_A_MATRIX, n, n, a,
IMSL_B_MATRIX, n, n, p_inva,
0);
/* Form a*g2_inverse(a)*a */
p_a_inva_a = imsl_f_mat_mul_rect("A*B",
IMSL_A_MATRIX, n, n, p_a_inva,
IMSL_B_MATRIX, n, n, a,
0);
imsl_f_write_matrix("The g2 inverse of a", n, n, p_inva,
0);
imsl_f_write_matrix("a*g2_inverse(a)\nviolates condition 3 of"
" the M-P inverse", n, n, p_a_inva,
0);
imsl_f_write_matrix("a = a*g2_inverse(a)*a\ncondition 1 of"
" the M-P inverse", n, n, p_a_inva_a,
0);
}
Output
The g2 inverse of a
1 2 3 4
1 0.0347 -0.0208 0.0000 0.0000
2 -0.0208 0.0903 0.0000 -0.0556
3 0.0000 0.0000 0.0000 0.0000
4 0.0000 -0.0556 0.0000 0.1111
a*g2_inverse(a)
violates condition 3 of the M-P inverse
1 2 3 4
1 1.0 -0.0 0.0 0.0
2 0.0 1.0 0.0 0.0
3 1.0 -0.5 0.0 0.0
4 0.0 -0.0 0.0 1.0
a = a*g2_inverse(a)*a
condition 1 of the M-P inverse
1 2 3 4
1 36 12 30 6
2 12 20 2 10
3 30 2 29 1
4 6 10 1 14
Warning Errors
IMSL_INCONSISTENT_EQUATIONS_2 | The linear system of equations is inconsistent.. |
IMSL_NOT_NONNEG_DEFINITE | The matrix A is not nonnegative definite. |
