Solves a system of linear equations Ax = b, with sparse complex coefficient matrix A. Using optional arguments, any of several related computations can be performed. These extra tasks include returning the LU factorization of A, computing the solution of Ax = b given an LU factorization, setting drop tolerances, and controlling iterative refinement.
Synopsis
#include <imsl.h>
f_complex *imsl_c_lin_sol_gen_coordinate (int n, int nz, Imsl_c_sparse_elem *a, f_complex *b, ..., 0)
The type double function is imsl_z_lin_sol_gen_coordinate.
Required Arguments
int n
(Input)
Number of rows in the matrix.
int nz
(Input)
Number of nonzeros in the matrix.
Imsl_c_sparse_elem *a
(Input)
Vector of length nz containing the
location and value of each nonzero entry in the matrix.
f_complex *b
(Input)
Vector of length n containing the
right-hand side.
Return Value
A pointer to the solution x of the sparse 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_coordinate (int n, int nz, Imsl_c_sparse_elem *a, f_complex *b,
IMSL_RETURN_SPARSE_LU_FACTOR, Imsl_c_sparse_lu_factor *lu_factor,
IMSL_SUPPLY_SPARSE_LU_FACTOR, Imsl_c_sparse_lu_factor *lu_factor,
IMSL_FREE_SPARSE_LU_FACTOR,
IMSL_RETURN_SPARSE_LU_IN_COORD, Imsl_c_sparse_elem **lu_coordinate, int **row_pivots, int **col_pivots,
IMSL_SUPPLY_SPARSE_LU_IN_COORD, int nzlu, Imsl_c_sparse_elem *lu_coordinate, int *row_pivots, int *col_pivots,
IMSL_FACTOR_ONLY,
IMSL_SOLVE_ONLY,
IMSL_RETURN_USER, f_complex x[],
IMSL_TRANSPOSE,
IMSL_CONDITION, float *condition,
IMSL_PIVOTING_STRATEGY, Imsl_pivot method,
IMSL_NUMBER_OF_SEARCH_ROWS, int num_search_row,
IMSL_ITERATIVE_REFINEMENT,
IMSL_DROP_TOLERANCE, float tolerance,
IMSL_HYBRID_FACTORIZATION, float density, int order_bound,
IMSL_GROWTH_FACTOR_LIMIT, float gf_limit,
IMSL_GROWTH_FACTOR, float *gf,
IMSL_SMALLEST_PIVOT, float *small_pivot
IMSL_NUM_NONZEROS_IN_FACTOR, int *num_nonzeros,
IMSL_CSC_FORMAT, int *col_ptr, int *row_ind, f_complex *values,
IMSL_MEMORY_BLOCK_SIZE, int block_size,
0)
Optional Arguments
IMSL_RETURN_SPARSE_LU_FACTOR,
Imsl_c_sparse_lu_factor *lu_factor
(Output)
The address of a structure of type Imsl_c_sparse_lu_factor.
The pointers within the structure are initialized to point to the LU
factorization by imsl_c_lin_sol_gen_coordinate.
IMSL_SUPPLY_SPARSE_LU_FACTOR,
Imsl_c_sparse_lu_factor *lu_factor
(Input)
The address of a structure of type Imsl_c_sparse_lu_factor.
This structure contains the LU factorization of the input matrix computed
by imsl_c_lin_sol_gen_coordinate
with the IMSL_RETURN_SPARSE_LU_FACTOR
option.
IMSL_FREE_SPARSE_LU_FACTOR,
Before
returning, free the linked list data structure containing the LU
factorization of A. Use this option only if the factors are no longer
required.
IMSL_RETURN_SPARSE_LU_IN_COORD,
Imsl_c_sparse_elem **lu_coordinate, int **row_pivots, int **col_pivots
(Output)
The LU factorization is returned in coordinate form in an
array of length nz in lu_coordinate. This is
more compact than the internal representation encapsulated in Imsl_c_sparse_lu. The
disadvantage is that during a SOLVE_ONLY call, the
internal representation of the factor must be reconstructed. If however, the
factor is to be stored after the program exits, and loaded again at some
subsequent run, the combination of IMSL_RETURN_LU_IN_COORD
and IMSL_SUPPLY_LU_IN_COORD
is probably the best choice, since the factors are in a format that is simple to
store and read.
IMSL_SUPPLY_SPARSE_LU_IN_COORD,
int nzlu,
Imsl_c_sparse_elem *lu_coordinate, int *row_pivots, int *col_pivots
(Input)
Supply the LU factorization in coordinate form. See IMSL_RETURN_SPARSE_LU_IN_COORD
for a description.
IMSL_FACTOR_ONLY,
Compute
the LU factorization of the input matrix and return. The argument b is ignored.
IMSL_SOLVE_ONLY,
Solve
Ax = b given the LU factorization of A.
This option requires the use of option IMSL_SUPPLY_SPARSE_LU_FACTOR
or IMSL_SUPPLY_SPARSE_LU_IN_COORD.
IMSL_RETURN_USER, f_complex x[]
(Output)
A user-allocated array of length n containing the solution
x.
IMSL_TRANSPOSE,
Solve
the problem ATx = b. This option
can be used in conjunction with either of the options that supply the
factorization.
IMSL_CONDITION, float *condition,
Estimate
the L1 condition number of A and return
in the variable condition.
IMSL_PIVOTING_STRATEGY,
Imsl_pivot method (Input)
Select the pivoting
strategy by setting method to one of the following: IMSL_ROW_MARKOWITZ,
IMSL_COLUMN_MARKOWITZ,
or IMSL_SYMMETRIC_MARKOWITZ.
Default:
IMSL_SYMMETRIC_MARKOWITZ.
IMSL_NUMBER_OF_SEARCH_ROWS, int num_search_row
(Input)
The number of rows which have the least number of nonzero elements
that will be searched for a pivot element.
Default: num_search_row = 3
IMSL_ITERATIVE_REFINEMENT,
Select
this option if iterative refinement is desired.
IMSL_DROP_TOLERANCE, float tolerance
(Input)
Possible fill-in is checked against tolerance. If the absolute value
of the new element is less than tolerance, it will be discarded.
Default:
tolerance = 0.0
IMSL_HYBRID_FACTORIZATION, float density, int order_bound,
(Input)
Enable the code to switch to a dense factorization method when the
density of the active submatrix reaches 0.0 ≤ density ≤ 1.0
and the order of the active submatrix is less than or equal to order_bound.
IMSL_GROWTH_FACTOR_LIMIT, float gf_limit
(Input)
The computation stops if the growth factor exceeds gf_limit.
Default: gf_limit = 1.e16
IMSL_GROWTH_FACTOR, float *gf
(Output)
gf
is calculated as the largest element in absolute value at any stage of the
Gaussian elimination divided by the largest element in absolute value in
A.
IMSL_SMALLEST_PIVOT, float *small_pivot
(Output)
A pointer to the value of the pivot element of smallest
magnitude.
IMSL_NUM_NONZEROS_IN_FACTOR, int *num_nonzeros
(Output)
A pointer to a scalar containing the total number of nonzeros in the
factor.
IMSL_CSC_FORMAT, int *col_ptr, int *row_ind, f_complex *values
(Input)
Accept the coefficient matrix in compressed sparse column (CSC)
format. See the main “Introduction” chapter at the beginning of this manual for
a discussion of this storage scheme.
IMSL_MEMORY_BLOCK_SIZE, int blocksize
(Input)
If space must be allocated for fill-in, allocate enough space for
blocksize new
nonzero elements.
Default: blocksize = nz
Description
The function imsl_c_lin_sol_gen_coordinate solves a system of linear equations Ax = b, where A is sparse. In its default use, it solves the so-called one off problem, by first performing an LU factorization of A using the improved generalized symmetric Markowitz pivoting scheme. The factor L is not stored explicitly because the saxpy operations performed during the elimination are extended to the right-hand side, along with any row interchanges. Thus, the system Ly = b is solved implicitly. The factor U is then passed to a triangular solver which computes the solution x from Ux = y.
If a sequence of systems Ax = b are to be solved where A is unchanged, it is usually more efficient to compute the factorization once, and perform multiple forward and back solves with the various right-hand sides. In this case the factor L is explicitly stored and a record of all row as well as column interchanges is made. The solve step then solves the two triangular systems Ly = b and Ux = y. The user specifies either the IMSL_RETURN_SPARSE_LU_FACTOR or the IMSL_RETURN_LU_IN_COORD option to retrieve the factorization, then calls the function subsequently with different right-hand sides, passing the factorization back in using either IMSL_SUPPLY_SPARSE_LU_FACTOR or IMSL_SUPPLY_SPARSE_LU_IN_COORD in conjunction with IMSL_SOLVE_ONLY. If IMSL_RETURN_SPARSE_LU_FACTOR is used, the final call to imsl_lin_sol_gen_coordinate should include IMSL_FREE_SPARSE_LU_FACTOR to release the heap used to store L and U.
If the solution to ATx = b is required, specify the option IMSL_TRANSPOSE. This keyword only alters the forward elimination and back substitution so that the operations UTy = b and LTx = y are performed to obtain the solution. So, with one call to produce the factorization, solutions to both Ax = b and ATx = b can be obtained.
The option IMSL_CONDITION is used to calculate and return an estimation of the L1 condition number of A. The algorithm used is due to Higham. Specification of IMSL_CONDITION causes a complete L to be computed and stored, even if a one off problem is being solved. This is due to the fact that Higham’s method requires solution to problems of the form Az = r and ATz = r.
The default pivoting strategy is symmetric Markowitz. If a row or column oriented problem is encountered, there may be some reduction in fill-in by selecting either IMSL_ROW_MARKOWITZ or IMSL_COLUMN_MARKOWITZ. The Markowitz strategy will search a pre-elected number of row or columns for pivot candidates. The default number is three, but this can be changed by using IMSL_NUM_OF_SEARCH_ROWS.
The option IMSL_DROP_TOLERANCE can be used to set a tolerance which can reduce fill-in. This works by preventing any new fill element which has magnitude less than the specified drop tolerance from being added to the factorization. Since this can introduce substantial error into the factorization, it is recommended that IMSL_ITERATIVE_REFINEMENT be used to recover more accuracy in the final solution. The trade-off is between space savings from the drop tolerance and the extra time needed in repeated solve steps needed for refinement.
The function imsl_c_lin_sol_gen_coordinate provides the option of switching to a dense factorization method at some point during the decomposition. This option is enabled by choosing IMSL_HYBRID_FACTORIZATION. One of the two parameters required by this option, density, specifies a minimum density for the active submatrix before a format switch will occur. A density of 1.0 indicates complete fill-in. The other parameter, order_bound, places an upper bound of the order of the active submatrix which will be converted to dense format. This is used to prevent a switch from occurring too early, possibly when the O(n3) nature of the dense factorization will cause performance degradation. Note that this option can significantly increase heap storage requirements.
Examples
Example 1
As an example, consider the following matrix:
Let
xT = (1 + i, 2 + 2i, 3 + 3i, 4 + 4i, 5 + 5i, 6 + 6i)
so that
Ax = (3 + 17i, −19 + 5i, 6 + 18i, − 38 + 32i, −63 + 49i, −57 + 83i)T
#include <imsl.h>
#include <stdlib.h>
int main()
{
static Imsl_c_sparse_elem a[] = {0, 0, {10.0, 7.0},
1, 1, {3.0, 2.0},
1, 2, {-3.0, 0.0},
1, 3, {-1.0, 2.0},
2, 2, {4.0, 2.0},
3, 0, {-2.0, -4.0},
3, 3, {1.0, 6.0},
3, 4, {-1.0, 3.0},
4, 0, {-5.0, 4.0},
4, 3, {-5.0, 0.0},
4, 4, {12.0, 2.0},
4, 5, {-7.0, 7.0},
5, 0, {-1.0, 12.0},
5, 1, {-2.0, 8.0},
5, 5, {3.0, 7.0}};
static f_complex b[] = {{3.0, 17.0}, {-19.0, 5.0}, {6.0, 18.0},
{-38.0, 32.0}, {-63.0, 49.0}, {-57.0, 83.0}};
int n = 6;
int nz = 15;
f_complex *x;
x = imsl_c_lin_sol_gen_coordinate (n, nz, a, b, 0);
imsl_c_write_matrix ("solution", n, 1, x, 0);
imsl_free (x);
}
Output
solution
1 ( 1, 1)
2 ( 2, 2)
3 ( 3, 3)
4 ( 4, 4)
5 ( 5, 5)
6 ( 6, 6)
Example 2
This examples sets A = E (1000, 10). A linear system is solved and the LU factorization returned. Then a second linear system is solved using the same coefficient matrix A just factored. Maximum absolute errors and execution time ratios are printed showing that forward and back solves take a small percentage of the computation time of a factor and solve. This ratio can vary greatly, depending on the order of the coefficient matrix, the initial number of nonzeros, and especially on the amount of fill-in produced during the elimination. Be aware that timing results are highly machine dependent.
#include <imsl.h>
#include <stdlib.h>
int main()
{
Imsl_c_sparse_elem *a;
Imsl_c_sparse_lu_factor lu_factor;
f_complex *b;
f_complex *x;
f_complex *mod_five;
f_complex *mod_ten;
float error_factor_solve;
float error_solve;
double time_factor_solve;
double time_solve;
int n = 1000;
int c = 10;
int i;
int nz;
int index;
/* Get the coefficient matrix */
a = imsl_c_generate_test_coordinate (n, c, &nz, 0);
/* Set two different predetermined solutions */
mod_five = (f_complex*) malloc (n*sizeof(*mod_five));
mod_ten = (f_complex*) malloc (n*sizeof(*mod_ten));
for (i=0; i<n; i++) {
mod_five[i] = imsl_cf_convert ((float)(i % 5), 0.0);
mod_ten[i] = imsl_cf_convert ((float)(i % 10), 0.0);
}
/* Choose b so that x will approximate mod_five */
b = imsl_c_mat_mul_rect_coordinate ("A*x",
IMSL_A_MATRIX, n, n, nz, a,
IMSL_X_VECTOR, n, mod_five,
0);
/* Time the factor/solve */
time_factor_solve = imsl_ctime();
x = imsl_c_lin_sol_gen_coordinate (n, nz, a, b,
IMSL_RETURN_SPARSE_LU_FACTOR, &lu_factor,
0);
time_factor_solve = imsl_ctime() - time_factor_solve;
/* Compute max abolute error */
error_factor_solve = imsl_c_vector_norm (n, x,
IMSL_SECOND_VECTOR, mod_five,
IMSL_INF_NORM, &index,
0);
imsl_free (b);
imsl_free (x);
/* Get new right hand side -- b = A * mod_ten */
b = imsl_c_mat_mul_rect_coordinate ("A*x",
IMSL_A_MATRIX, n, n, nz, a,
IMSL_X_VECTOR, n, mod_ten,
0);
/* Use the previously computed factorization
to solve Ax = b */
time_solve = imsl_ctime();
x = imsl_c_lin_sol_gen_coordinate (n, nz, a, b,
IMSL_SUPPLY_SPARSE_LU_FACTOR, &lu_factor,
IMSL_SOLVE_ONLY,
0);
time_solve = imsl_ctime() - time_solve;
error_solve = imsl_c_vector_norm (n, x,
IMSL_SECOND_VECTOR, mod_ten,
IMSL_INF_NORM, &index,
0);
imsl_free (b);
imsl_free (x);
/* Print errors and ratio of execution times */
printf ("absolute error (factor/solve) = %e\n",
error_factor_solve);
printf ("absolute error (solve) = %e\n", error_solve);
printf ("time_solve/time_factor_solve = %f\n",
time_solve/time_factor_solve);
}
Output
absolute error (factor/solve) = 2.389053e-06
absolute error (solve) = 7.656095e-06
time_solve/time_factor_solve = 0.070313
