Solves a real 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 ATx = b, or computing the solution of Ax = b given the LU factorization of A.
Synopsis
#include <imsl.h>
float *imsl_f_lin_sol_gen (int n, float a[], float b[], …, 0)
The type double procedure is imsl_d_lin_sol_gen.
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 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_gen (int
n, float
a[],
float
b[],
IMSL_A_COL_DIM, int
a_col_dim,
IMSL_TRANSPOSE,
IMSL_RETURN_USER, float x[],
IMSL_FACTOR, int
**p_pvt, float
**p_factor,
IMSL_FACTOR_USER, int
pvt[],
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
inva_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_TRANSPOSE
Solve
ATx = b.
Default:
Solve
Ax = b
IMSL_RETURN_USER, float x[]
(Output)
A user-allocated array of length n containing the solution
x.
IMSL_FACTOR, int **p_pvt, float **p_factor (Output)
p_pvt: 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_f_lin_sol_gen. Typically, int *p_pvt is declared, and &p_pvt is used as an argument.
p_factor: 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_f_lin_sol_gen. The lower-triangular part of this array contains information necessary to construct L, and the upper-triangular part contains U. Typically, float *p_factor is declared, and &p_factor is used as an argument.
IMSL_FACTOR_USER, int pvt[], float factor[] (Input/Output)
pvt[]: A user-allocated array of size n containing the pivot sequence for the factorization.
factor[]: A user-allocated array of size n × n containing the LU factorization of A. The strictly lower-triangular part of this array contains information necessary to construct L, and the upper-triangular part contains U. If A is not needed, factor and a can share the same storage.
These parameters are input if IMSL_SOLVE is specified. They are output otherwise.
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, 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_gen.
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. This option cannot be used with the option IMSL_SOLVE_ONLY.
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_f_lin_sol_gen is
NULL.
IMSL_SOLVE_ONLY
Solve
Ax = b given the LU factorization previously
computed by imsl_f_lin_sol_gen. By
default, the solution to Ax = b is pointed to by imsl_f_lin_sol_gen. 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_gen is
NULL.
Description
The function imsl_f_lin_sol_gen
solves a system of linear algebraic equations with a real coefficient matrix
A. It first computes the LU factorization of A with partial
pivoting such that
L-1A = U.
The matrix U is upper triangular, while L-1A ≡ Pn Ln-nPn-1… L1P1A ≡ U. The factors
Pi and Li are defined by the
partial pivoting. Each Pi is an interchange of
row i with row j ³ i. Thus, Pi is defined by that
value of j. Every
is an elementary elimination matrix. The vector mi is zero in entries 1, …, i. This vector is stored as column i in the strictly lower-triangular part of the working array containing the decomposition information.
The factorization efficiency is based on a technique of “loop unrolling and jamming” by Dr. Leonard J. Harding of the University of Michigan, Ann Arbor, Michigan. The solution of the linear system is then found by solving two simpler systems, y = L-1b and x = U-1y. 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_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. This is the simplest use of the function. The equations follow below:
x1 + 3x2 + 3x3 = 1
x1 + 3x2 + 4x3 = 4
x1 + 4x2 + 3x3 = −1
#include
<imsl.h>
main()
{
int n = 3;
float *x;
float a[] = {1.0, 3.0,
3.0,
1.0, 3.0,
4.0,
1.0, 4.0, 3.0};
float
b[] = {1.0, 4.0,
-1.0};
/* Solve Ax = b
for x */
x = imsl_f_lin_sol_gen (n, a, b,
0);
/* Print x */
imsl_f_write_matrix ("Solution, x, of Ax =
b", 1, 3, x, 0);
}
Output
Solution, x, of Ax =
b
1
2
3
-2
-2
3
Example 2
This example solves the transpose problem ATx = b and returns the LU factorization of A with partial pivoting. The same data as the initial example is used, except the solution x = A-Tb is returned in an array allocated in the main program. The L matrix is returned in implicit form.
#include
<imsl.h>
main()
{
int n = 3,
pvt[3];
float
factor[9];
float
x[3];
float a[]
= {1.0, 3.0,
3.0,
1.0, 3.0,
4.0,
1.0,
4.0, 3.0};
float
b[] = {1.0, 4.0,
-1.0};
/* Solve trans(A)*x = b for x */
imsl_f_lin_sol_gen (n, a,
b,
IMSL_TRANSPOSE,
IMSL_RETURN_USER,
x,
IMSL_FACTOR_USER, pvt,
factor,
0);
/* Print x */
imsl_f_write_matrix ("Solution, x, of
trans(A)x = b", 1, n, x,
0);
/* Print factors and pivot sequence */
imsl_f_write_matrix
("LU factors of A", n, n, factor, 0);
imsl_i_write_matrix
("Pivot sequence", 1, n, pvt, 0);
}
Output
Solution, x, of trans(A)x = b
1
2
3
4
-4
1
LU
factors of
A
1
2
3
1
1
3
3
2
-1
1
0
3
-1
0
1
Pivot sequence
1 2
3
1 3 3
Example 3
This example computes the inverse of the 3 × 3 matrix A of the initial example and solves the same linear system. The matrix product C = A-1A is computed and printed. The function imsl_f_mat_mul_rect is used to compute C. The approximate result C = I is obtained.
#include <imsl.h>
float
a[] = {1.0, 3.0,
3.0,
1.0, 3.0,
4.0,
1.0, 4.0, 3.0};
float b[] = {1.0, 4.0,
-1.0};
main()
{
int n =
3;
float
*x;
float
*p_inva;
float
*C;
/* Solve Ax = b
*/
x = imsl_f_lin_sol_gen (n, a,
b,
IMSL_INVERSE,
&p_inva,
0);
/* Print solution */
imsl_f_write_matrix ("Solution,
x, of Ax = b", 1, n, x,
0);
/* Print input and
inverse matrices */
imsl_f_write_matrix ("Input A", n, n,
a, 0);
imsl_f_write_matrix ("Inverse of A", n, n, p_inva,
0);
/* Check result and print */
C =
imsl_f_mat_mul_rect("A*B",
IMSL_A_MATRIX, n, n, p_inva,
IMSL_B_MATRIX, n, n, a,
0);
imsl_f_write_matrix ("Product matrix,
inv(A)*A",n,n,C,0);
}
Output
Solution, x, of Ax = b
1
2 3
-2
-2
3
Input A
1
2
3
1
1
3
3
2
1
3
4
3
1
4
3
Inverse of
A
1
2
3
1
7
-3
-3
2
-1
0
1
3
-1
1
0
Product matrix,
inv(A)*A
1
2
3
1
1
0
0
2
0
1
0
3
0
0 1
Example 4
This example computes the solution of two systems. Only the right-hand sides differ. The matrix and first right-hand side are given in the initial example. The second right-hand side is the vector c = [0.5, 0.3, 0.4]T. The factorization information is computed with the first solution and is used to compute the second solution. The factorization work done in the first step is avoided in computing the second solution.
#include
<imsl.h>
main()
{
int n = 3,
pvt[3];
float
factor[9];
float
*x,*y;
float a[]
= {1.0, 3.0,
3.0,
1.0, 3.0,
4.0,
1.0, 4.0, 3.0};
float b[] = {1.0, 4.0,
-1.0};
float c[]
= {0.5, 0.3,
0.4};
/*
Solve A*x = b for x */
x =
imsl_f_lin_sol_gen (n, a,
b,
IMSL_FACTOR_USER, pvt,
factor,
0);
/* Print x */
imsl_f_write_matrix ("Solution, x, of Ax =
b", 1, n, x,
0);
/* Solve for A*y = c for y */
y = imsl_f_lin_sol_gen (n,
a,
c,
IMSL_SOLVE_ONLY,
IMSL_FACTOR_USER, pvt,
factor,
0);
imsl_f_write_matrix ("Solution, y, of Ay = c", 1, n,
y, 0);
}
Output
Solution, x, of Ax = b
1
2 3
-2
-2
3
Solution, y, of Ay = c
1
2
3
1.4
-0.1 -0.2
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_SINGULAR_MATRIX The input matrix is singular.
|
Visual Numerics, Inc. PHONE: 713.784.3131 FAX:713.781.9260 |