lin_sol_posdef (complex)
Solves a complex Hermitian 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 = LLH or computing the solution to Ax = b given the Cholesky factor, L.
Synopsis
#include <imsl.h>
f_complex *imsl_c_lin_sol_posdef (int n, f_complex a[], f_complex b[], …, 0)
The type d_complex function is imsl_z_lin_sol_posdef.
Required Arguments
int n (Input)
Number of rows and columns in the matrix.
f_complex a[] (Input)
Array of size n × n containing 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 Hermitian positive definite 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 (int n, f_complex a[], f_complex b[],
IMSL_A_COL_DIM, int a_col_dim,
IMSL_RETURN_USER, f_complex x[],
IMSL_FACTOR, f_complex **p_factor,
IMSL_FACTOR_USER, f_complex factor[],
IMSL_FAC_COL_DIM, int fac_col_dim,
IMSL_CONDITION, float *cond,
IMSL_FACTOR_ONLY,
IMSL_SOLVE_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, f_complex x[] (Output)
A user-allocated array of size n containing the solution x.
IMSL_FACTOR, f_complex **p_factor (Output)
The address of a pointer to an array of size n × n containing the LLH factorization of A. On return, the necessary space is allocated by imsl_c_lin_sol_posdef. The lower‑triangular part of this array contains L, and the upper-triangular part contains LH. 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 n × n containing the LLH factorization of A. The lower‑triangular part of this array contains L, and the upper-triangular part contains LH. If A is not needed, a and factor can share the same storage. If IMSL_SOLVE is specified, factor is input. Otherwise, it is output.
IMSL_FAC_COL_DIM, int fac_col_dim (Input)
The column dimension of the array containing the LLH factorization of A.
Default: fac_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 LLH 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 is NULL.
IMSL_SOLVE_ONLY
Solve Ax = b given the LLH factorization previously computed by imsl_c_lin_sol_posdef. By default, the solution to Ax = b is pointed to by imsl_c_lin_sol_posdef. If IMSL_SOLVE_ONLY is used, argument IMSL_FACTOR_USER is required and argument a is ignored.
Description
The function imsl_c_lin_sol_posdef solves a system of linear algebraic equations having a Hermitian positive definite coefficient matrix A. The function first computes the LLH factorization of A. The solution of the linear system is then found by solving the two simpler systems, y = L-1b and x = L-Hy. When the solution to the linear system is required, an estimate of the L1 condition number of A is computed using the algorithm 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_c_lin_sol_posdef fails if L, the lower-triangular matrix in the factorization, has a zero diagonal element.
Examples
Example 1
A system of five linear equations with a Hermitian positive definite coefficient matrix is solved in this example. The equations are as follows:
2x1 +(−1 + i)x2 = 1 +5i
(−1 − i)x1 +4x2 + (1 + 2i)x3 = 12 − 6i
(1 − 2i)x2 +10x3 + 4ix4 = 1 − 16i
−4ix3 + 6x4 + (1 + i)x5 = −3 − 3i
(1 − i)x4 + 9x5 = 25 + 16i
#include <imsl.h>
int main()
{
int n = 5;
f_complex *x;
f_complex a[] = {
{2.0,0.0}, {-1.0,1.0},{0.0,0.0}, {0.0,0.0}, {0.0,0.0},
{-1.0,-1.0},{4.0,0.0}, {1.0,2.0}, {0.0,0.0}, {0.0,0.0},
{0.0,0.0}, {1.0,-2.0},{10.0,0.0},{0.0,4.0}, {0.0,0.0},
{0.0,0.0}, {0.0,0.0}, {0.0,-4.0},{6.0,0.0}, {1.0,1.0},
{0.0,0.0}, {0.0,0.0}, {0.0,0.0}, {1.0,-1.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}
};
/* Solve Ax = b for x */
x = imsl_c_lin_sol_posdef(n, a, b, 0);
/* Print x */
imsl_c_write_matrix("Solution, x, of Ax = b", 1, n, x, 0);
}
Output
Solution, x, of Ax = b
1 2 3
( 2, 1) ( 3, -0) ( -1, -1)
4 5
( 0, -2) ( 3, 2)
Example 2
This example solves the same system of five linear equations as in the first example. This time, the LLH factorization of A and the solution x is returned in an array allocated in the main program.
#include <imsl.h>
int main()
{
int n = 5;
f_complex x[5], *p_factor;
f_complex a[] = {
{2.0,0.0}, {-1.0,1.0},{0.0,0.0}, {0.0,0.0}, {0.0,0.0},
{-1.0,-1.0},{4.0,0.0}, {1.0,2.0}, {0.0,0.0}, {0.0,0.0},
{0.0,0.0}, {1.0,-2.0},{10.0,0.0},{0.0,4.0}, {0.0,0.0},
{0.0,0.0}, {0.0,0.0}, {0.0,-4.0},{6.0,0.0}, {1.0,1.0},
{0.0,0.0}, {0.0,0.0}, {0.0,0.0}, {1.0,-1.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}
};
/* Solve Ax = b for x */
imsl_c_lin_sol_posdef(n, a, b,
IMSL_RETURN_USER, x,
IMSL_FACTOR, &p_factor,
0);
/* Print x */
imsl_c_write_matrix("Solution, x, of Ax = b", 1, n, x, 0);
/* Print Cholesky factor of A */
imsl_c_write_matrix("Cholesky factor L, and ctrans(L), of A",
n, n, p_factor, 0);
}
Output
Solution, x, of Ax = b
1 2 3
( 2, 1) ( 3, -0) ( -1, -1)
4 5
( 0, -2) ( 3, 2)
Cholesky factor L, and ctrans(L), of A
1 2 3
1 ( 1.414, 0.000) ( -0.707, 0.707) ( 0.000, -0.000)
2 ( -0.707, -0.707) ( 1.732, 0.000) ( 0.577, 1.155)
3 ( 0.000, 0.000) ( 0.577, -1.155) ( 2.887, 0.000)
4 ( 0.000, 0.000) ( 0.000, 0.000) ( 0.000, -1.386)
5 ( 0.000, 0.000) ( 0.000, 0.000) ( 0.000, 0.000)
4 5
1 ( 0.000, -0.000) ( 0.000, -0.000)
2 ( 0.000, -0.000) ( 0.000, -0.000)
3 ( 0.000, 1.386) ( 0.000, -0.000)
4 ( 2.020, 0.000) ( 0.495, 0.495)
5 ( 0.495, -0.495) ( 2.917, 0.000)
Warning Errors
IMSL_HERMITIAN_DIAG_REAL_1 | The diagonal of a Hermitian matrix must be real. Its imaginary part is set to zero. |
IMSL_HERMITIAN_DIAG_REAL_2 | The diagonal of a Hermitian matrix must be real. The imaginary part will be used as zero in the algorithm. |
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 # minor matrix of the input matrix is not positive definite. |
IMSL_HERMITIAN_DIAG_REAL | During the factorization the matrix has a large imaginary component on the diagonal. Thus, it cannot be positive definite. |
IMSL_SINGULAR_TRI_MATRIX | The triangular matrix is singular. The index of the first zero diagonal term is #. |
