lin_sol_gen (complex)

CNLMath : Linear Systems : lin_sol_gen (complex)

Synopsis

Required Arguments

Return Value

Synopsis with Optional Arguments

more...

Solves a complex general system of linear equations Ax = b. Using optional arguments, any of several related computations can be performed. These extra tasks include computing the LU factorization of A using partial pivoting, computing the inverse matrix A-1, solving AHx = b, or computing the solution of Ax = b given the LU factorization of A.

Synopsis

#include <imsl.h>

f_complex *imsl_c_lin_sol_gen (int n, f_complex a[], f_complex b[], …, 0)

The type d_complex function is imsl_z_lin_sol_gen.

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

Synopsis with Optional Arguments

#include <imsl.h>

f_complex *imsl_c_lin_sol_gen (int n, f_complex a[], f_complex b[],

IMSL_A_COL_DIM, int a_col_dim,

IMSL_TRANSPOSE,

IMSL_RETURN_USER, f_complex x[],

IMSL_FACTOR, int **p_pvt, f_complex **p_factor,

IMSL_FACTOR_USER, int pvt[], f_complex factor[],

IMSL_FAC_COL_DIM, int fac_col_dim,

IMSL_INVERSE, f_complex **p_inva,

IMSL_INVERSE_USER, f_complex inva[],

IMSL_INV_COL_DIM, int inva_col_dim,

IMSL_CONDITION, float *cond,

IMSL_ITERATIVE_REFINEMENT, int refine,

IMSL_FACTOR_ONLY,

IMSL_SOLVE_ONLY,

IMSL_INVERSE_ONLY,

Optional Arguments

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

IMSL_TRANSPOSE
Solve AHx = b
Default: Solve Ax = b

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

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

int **p_pvt (Output)
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_c_lin_sol_gen. Typically, int *p_pvt is declared, and &p_pvt is used as an argument.

f_complex **p_factor (Output)
The address of a pointer to an array of size n × n containing the LU factorization of A with column pivoting. On return, the necessary space is allocated by imsl_c_lin_sol_gen. The lower-triangular part of this array contains information necessary to construct L, and the upper-triangular part contains U. Typically, f_complex *p_factor is declared, and &p_factor is used as an argument.

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

int pvt[] (Input/Output)
A user-allocated array of size n containing the pivot sequence for the factorization.

f_complex factor[] (Input/Output)
A user-allocated array of size n × n containing the LU factorization of A. The lower-triangular part of this array contains information necessary to construct L, and the upper-triangular part contains U.

These parameters are input if IMSL_SOLVE is specified. They are output otherwise. If A is not needed, factor and a can share the same storage.

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

IMSL_INVERSE, f_complex **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_c_lin_sol_gen. Typically, f_complex *p_inva is declared, and &p_inva is used as an argument.

IMSL_INVERSE_USER, f_complex 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_ITERATIVE_REFINEMENT, int refine (Input)
Indicates if iterative refinement is desired.

refine	Description
0	No iterative refinement.
1	Do iterative refinement.

Default: refine= 0.

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_c_lin_sol_gen is NULL.

IMSL_SOLVE_ONLY
Solve Ax = b given the LU factorization previously computed by imsl_c_lin_sol_gen. By default, the solution to Ax = b is pointed to by imsl_c_lin_sol_gen. If IMSL_SOLVE_ONLY is used, argument IMSL_FACTOR_USER is required. If iterative refinement of the solution is desired, argument a must be present. Otherwise, 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. Argument b is then ignored, and the returned value of imsl_c_lin_sol_gen is NULL.

Description

The function imsl_c_lin_sol_gen solves a system of linear algebraic equations with a complex coefficient matrix A. It first computes the LU factorization of A with partial pivoting such that L-1A = U. Let F be the matrix p_factor returned by optional argument IMSL_FACTOR. The triangular matrix U is stored in the upper triangle of F. The strict lower triangle of F contains the information needed to reconstruct L‑1 using

The factors Pi and Li are defined by partial pivoting. Pi is the identity matrix with rows i and p_pvt[i-1] interchanged. Li is the identity matrix with Fji , for j = i + 1,…, n, inserted below the diagonal in column i.

The solution of the linear system is then found by solving two simpler systems, y = L-1b and x = U -1y. Additionally, the accuracy of the solution can be improved by iterative refinement. IMSL uses mixed precision iterative refinement in single precision and fixed precision iterative refinement in double precision. In double precision, the residuals b-Ax are computed with high accuracy using algorithms based on Ogita, Rump and Oishi (2005). When the solution to the linear system or the inverse of the matrix is computed, 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_c_lin_sol_gen fails if U, the upper-triangular part of the factorization, has a zero diagonal element.

Examples

Example 1

This example solves a system of three linear equations. The equations are:

(1 + i) x1 + (2 + 3i) x2 + (3 − 3i) x3 = 3 + 5i

(2 + i) x1 + (5 + 3i) x2 + (7 − 5i) x3 = 22 + 10i

(−2 + i) x1 + (−4 + 4i) x2 + (5 + 3i) x3 = −10 + 4i

#include <imsl.h>

f_complex a[] = {{1.0, 1.0}, {2.0, 3.0}, {3.0, -3.0},

{2.0, 1.0}, {5.0, 3.0}, {7.0, -5.0},

{-2.0, 1.0}, {-4.0, 4.0}, {5.0, 3.0}};

f_complex b[] = {{3.0, 5.0}, {22.0, 10.0}, {-10.0, 4.0}};

int main()

{

int n = 3;

f_complex *x;

/* Solve Ax = b for x */

x = imsl_c_lin_sol_gen (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

( 1, -1) ( 2, 4) ( 3, -0)

Example 2

This example solves the conjugate transpose problem AHx = b and returns the LU factorization of A using partial pivoting. This example differs from the first example in that the solution array is allocated in the main program.

#include <imsl.h>

f_complex a[] = {{1.0, 1.0}, {2.0, 3.0}, {3.0, -3.0},

{2.0, 1.0}, {5.0, 3.0}, {7.0, -5.0},

{-2.0, 1.0}, {-4.0, 4.0}, {5.0, 3.0}};

f_complex b[] = {{3.0, 5.0}, {22.0, 10.0}, {-10.0, 4.0}};

int main()

{

int n = 3, pvt[3];

f_complex factor[9];

f_complex x[3];

/* Solve ctrans(A)*x = b for x */

imsl_c_lin_sol_gen (n, a, b,

IMSL_TRANSPOSE,

IMSL_RETURN_USER, x,

IMSL_FACTOR_USER, pvt, factor,

0);

/* Print x */

imsl_c_write_matrix ("Solution, x, of ctrans(A)x = b", 1, n, x, 0);

/* Print factors and pivot sequence */

imsl_c_write_matrix ("LU factors of A", n, n, factor, 0);

imsl_i_write_matrix ("Pivot sequence", 1, n, pvt, 0);

}

Output

Solution, x, of ctrans(A)x = b

1 2 3

( -9.79, 11.23) ( 2.96, -3.13) ( 1.85, 2.47)

LU factors of A

1 2 3

1 ( -2.000, 1.000) ( -4.000, 4.000) ( 5.000, 3.000)

2 ( 0.600, 0.800) ( -1.200, 1.400) ( 2.200, 0.600)

3 ( 0.200, 0.600) ( -1.118, 0.529) ( 4.824, 1.294)

Pivot sequence

1 2 3

3 3 3

Example 3

This example computes the inverse of the 3 × 3 matrix A in the first example and also solves the linear system. The product matrix C = A-1A is computed as a check. The approximate result is C = I.

#include <imsl.h>

f_complex a[] = {{1.0, 1.0}, {2.0, 3.0}, {3.0, -3.0},

{2.0, 1.0}, {5.0, 3.0}, {7.0, -5.0},

{-2.0, 1.0}, {-4.0, 4.0}, {5.0, 3.0}};

f_complex b[] = {{3.0, 5.0}, {22.0, 10.0}, {-10.0, 4.0}};

int main()

{

int n = 3;

f_complex *x;

f_complex *p_inva;

f_complex *C;

/* Solve Ax = b for x */

x = imsl_c_lin_sol_gen (n, a, b,

IMSL_INVERSE, &p_inva,

0);

/* Print solution */

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

/* Print input and inverse matrices */

imsl_c_write_matrix ("Input A", n, n, a, 0);

imsl_c_write_matrix ("Inverse of A", n, n, p_inva, 0);

/* Check and print result */

C = imsl_c_mat_mul_rect ("A*B",

IMSL_A_MATRIX, n,n, p_inva,

IMSL_B_MATRIX, n,n, a,

0);

imsl_c_write_matrix ("Product, inv(A)*A", n, n, C, 0);

}

Output

Solution, x, of Ax = b

1 2 3

( 1, -1) ( 2, 4) ( 3, -0)

Input A

1 2 3

1 ( 1, 1) ( 2, 3) ( 3, -3)

2 ( 2, 1) ( 5, 3) ( 7, -5)

3 ( -2, 1) ( -4, 4) ( 5, 3)

Inverse of A

1 2 3

1 ( 1.330, 0.594) ( -0.151, 0.028) ( -0.604, 0.613)

2 ( -0.632, -0.538) ( 0.160, 0.189) ( 0.142, -0.245)

3 ( -0.189, 0.160) ( 0.193, -0.052) ( 0.024, 0.042)

Product, inv(A)*A

1 2 3

1 ( 1, -0) ( -0, -0) ( -0, 0)

2 ( 0, 0) ( 1, 0) ( 0, -0)

3 ( -0, -0) ( -0, 0) ( 1, 0)

Warning Errors

IMSL_ILL_CONDITIONED	The input matrix is too ill-conditioned. An estimate of the reciprocal of the L1 condition number is “rcond” = #. The solution might not be accurate.
IMSL_ILL_CONDITIONED_1	The input matrix is too ill-conditioned for iterative refinement to be effective.

Fatal Errors

IMSL_SINGULAR_MATRIX

The input matrix is singular.