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 return­ing 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 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,
                              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_NUM_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. 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, Imsl_c_sparse_elem *lu_coordinate,
int *row_pivots, int *col_pivots   (Output)
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_NUM_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,
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_FACTOR_RESIZE_INCREMENT, int increment   (Input)
Supply the number of nonzeros which will be added to the factor if current allocations are inadequate.
Default: increment = 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 elimina­tion 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 effi­cient to compute the factorization once, and perform multiple forward and back solves with the var­ious 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_RE­TURN_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_S­PARSE_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 condi­tion number of A. The algorithm used is due to Higham. Specification of IMSL_CONDITION caus­es 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, by 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 factoriza­tion will cause performance degradation. Note that this option can significantly increase heap stor­age 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>

 

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);

 

        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 re­turned. 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 espe­cially on the amount of fill-in produced during the elimination. Be aware that timing results are high­ly machine dependent.

#include <imsl.h>

#include <stdlib.h>

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);

        free (b);

        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);

        free (b);

        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

 

 

Visual Numerics, Inc. Visual Numerics - Developers of IMSL and PV-WAVE http://www.vni.com/ PHONE: 713.784.3131 FAX:713.781.9260