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.
#include <imsl.h>
f_complex *imsl_c_lin_sol_posdef (int n, f_complex a[], f_complex b[], …, 0)
The type d_complex procedure is imsl_z_lin_sol_posdef.
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.
A pointer to the solution x of the Hermitian positive definite linear system Ax = b. To release this space, use free. If no solution was computed, then NULL is returned.
#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)
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.
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-1y. 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.
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>
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);
}
Solution, x, of Ax =
b
1
2
3
(
2, 1)
(
3, -0)
(
-1,
-1)
4
5
(
0, -2)
(
3, 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>
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);
}
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)
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.
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 #.
|
Visual Numerics, Inc. PHONE: 713.784.3131 FAX:713.781.9260 |