Solves a sparse real symmetric positive definite system of linear equations A = b. Using optional arguments, any of several related computations can be performed. These extra tasks include returning the symbolic factorization of A, returning the numeric factorization of A, and computing the solution of Ax = b given either the symbolic or numeric factorizations.
Synopsis
#include <imsl.h>
float *imsl_f_lin_sol_posdef_coordinate (int n, int nz, Imsl_f_sparse_elem *a, float *b, ..., 0)
The type double function is imsl_d_lin_sol_posdef_coordinate.
Required Arguments
int n
(Input)
Number of rows in the matrix.
int nz
(Input)
Number of nonzeros in lower triangle of the matrix.
Imsl_f_sparse_elem *a
(Input)
Vector of length nz containing the
location and value of each nonzero entry in the lower triangle of the
matrix.
float *b
(Input)
Vector of length n containing the
right-hand side.
Return Value
A pointer to the solution x of the sparse symmetric positive definite 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>
float *imsl_f_lin_sol_posdef_coordinate (int n, int nz, Imsl_f_sparse_elem *a, float *b,
IMSL_RETURN_SYMBOLIC_FACTOR, Imsl_symbolic_factor *sym_factor,
IMSL_SUPPLY_SYMBOLIC_FACTOR, Imsl_symbolic_factor *sym_factor,
IMSL_SYMBOLIC_FACTOR_ONLY,
IMSL_RETURN_NUMERIC_FACTOR, Imsl_f_numeric_factor *num_factor,
IMSL_SUPPLY_NUMERIC_FACTOR, Imsl_f_numeric_factor *num_factor,
IMSL_NUMERIC_FACTOR_ONLY,
IMSL_SOLVE_ONLY,
IMSL_MULTIFRONTAL_FACTORIZATION,
IMSL_RETURN_USER, float x[],
IMSL_SMALLEST_DIAGONAL_ELEMENT, float *small_element,
IMSL_LARGEST_DIAGONAL_ELEMENT, float *largest_element,
IMSL_NUM_NONZEROS_IN_FACTOR, int *num_nonzeros,
IMSL_CSC_FORMAT, int *col_ptr, int *row_ind, float *values,
0)
Optional Arguments
IMSL_RETURN_SYMBOLIC_FACTOR,
Imsl_symbolic_factor *sym_factor
(Output)
A pointer to a structure of type Imsl_symbolic_factor
containing, on return, the symbolic factorization of the input matrix.
IMSL_SUPPLY_SYMBOLIC_FACTOR,
Imsl_symbolic_factor *sym_factor
(Input)
A pointer to a structure of type Imsl_symbolic_factor. This
structure contains the symbolic factorization of the input matrix computed by
imsl_f_lin_sol_posdef_coordinate
with the IMSL_RETURN_SYMBOLIC_FACTOR
option.
IMSL_SYMBOLIC_FACTOR_ONLY,
Compute
the symbolic factorization of the input matrix and return. The argument b is ignored.
IMSL_RETURN_NUMERIC_FACTOR,
Imsl_f_numeric_factor *num_factor
(Output)
A pointer to a structure of type Imsl_f_numeric_factor
containing, on return, the numeric factorization of the input matrix. A detailed
description of the Imsl_f_numeric_factor structure is given in the
following table:
|
Parameter |
Data |
Description |
|
nzsub |
int ** |
A pointer to an array containing the row subscripts for the non-zero off-diagonal elements of the Cholesky factor. This array is allocated to be of length nz but all elements of the array may not be used. |
|
xnzsub |
int ** |
A pointer to an array of length n + 1 containing indices for nzsub. The row subscripts for the non-zeros in column j of the cholesky factor are stored consecutively beginning with nzsub[xnzsub[j]]. |
|
xlnz |
int ** |
A pointer to an array of length n + 1 containing the starting and stopping indices to use to extract the non-zero off-diagonal elements from array alnz. For column j of the factor matrix, the starting and stopping indices of alnz are stored in xlnz[j] and xlnz[j + 1] respectively. |
|
alnz |
float ** |
A pointer to an array of length nz - n containing the non-zero off-diagonal elements of the Cholesky factor. |
|
perm |
int ** |
A pointer to an array of length n containing the permutation vector. |
|
diag |
float ** |
A pointer to an array of length n containing the diagonal elements of the Cholesky factor. |
Let L be the numeric factorization of a. In the structure
described above, the diagonal elements of L are stored in diag. The off-diagonal
non-zero elements of L are stored in alnz. The starting and
stopping indices to use to extract the non-zero elements of L from alnz for column
j are stored in xlnz[j] and xlnz[j + 1] respectively. The
row indices of the non-zero elements of
L are contained in nzsub. xnzsub[i] contains the index
of nzsub from
which one should start to extract the row indices for L for column
i. This is best illustrated by the following code fragment which
reconstructs the lower triangle of the factor matrix L from the
components of the above structure:
Imsl_f_numeric_factor
numfctr;
.
.
.
for (i = 0; i < n; i++){
L[i][i] = (*numfctr.diag)[i];
if((*numfctr.xlnz)[i] >
(nz-n)+1) break;
start =
(*numfctr.xlnz)[i]-1;
stop =
(*numfctr.xlnz)[i+1]-1;
k =
(*numfctr.xnzsub)[i]-1;
for (j = start; j < stop;
j++){
L[(*numfctr.nzsub)[k]-1][i]
= (*numfctr.alnz)[j];
k++;
}
}
IMSL_SUPPLY_NUMERIC_FACTOR,
Imsl_f_numeric_factor *num_factor
(Input)
A pointer to a structure of type Imsl_f_numeric_factor. This
structure contains the numeric factorization of the input matrix computed by
imsl_f_lin_sol_posdef_coordinate
with the IMSL_RETURN_NUMERIC_FACTOR
option. The structure is described in the IMSL_RETURN_NUMERIC_FACTOR
optional argument desription.
IMSL_NUMERIC_FACTOR_ONLY,
Compute
the numeric factorization of the input matrix and return. The argument b is
ignored.
IMSL_SOLVE_ONLY,
Solve
Ax = b given the numeric or symbolic factorization of
A. This option requires the use of either IMSL_SUPPLY_NUMERIC_FACTOR
or IMSL_SUPPLY_SYMBOLIC_FACTOR.
IMSL_MULTIFRONTAL_FACTORIZATION,
Perform
the numeric factorization using a multifrontal technique. By default, a standard
factorization is computed based on a sparse compressed storage scheme.
IMSL_RETURN_USER, float x[]
(Output)
A user-allocated array of length n containing the solution
x.
IMSL_SMALLEST_DIAGONAL_ELEMENT, float *small_element
(Output)
A pointer to a scalar containing the smallest diagonal element that
occurred during the numeric factorization. This option is valid only if the
numeric factorization is computed during this call to imsl_f_lin_sol_posdef_coordinate.
IMSL_LARGEST_DIAGONAL_ELEMENT, float *large_element
(Output)
A pointer to a scalar containing the largest diagonal element that
occurred during the numeric factorization. This option is valid only if the
numeric factorization is computed during this call to imsl_f_lin_sol_posdef_coordinate.
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, float *values
(Input)
Accept the coefficient matrix in compressed sparse column (CSC)
format. See the “Matrix Storage Modes” section of the “Introduction” at the
beginning of this manual for a discussion of this storage scheme.
Description
The function imsl_f_lin_sol_posdef_coordinate solves a system of linear algebraic equations having a sparse symmetric positive definite coefficient matrix A. In this function’s default usage, a symbolic factorization of a permutation of the coefficient matrix is computed first. Then a numerical factorization is performed. The solution of the linear system is then found using the numeric factor.
The symbolic factorization step of the computation consists of determining a minimum degree ordering and then setting up a sparse data structure for the Cholesky factor, L. This step only requires the “pattern” of the sparse coefficient matrix, i.e., the locations of the nonzeros elements but not any of the elements themselves. Thus, the val field in the Imsl_f_sparse_elem structure is ignored. If an application generates different sparse symmetric positive definite coefficient matrices that all have the same sparsity pattern, then by using IMSL_RETURN_SYMBOLIC_FACTOR and IMSL_SUPPLY_SYMBOLIC_FACTOR, the symbolic factorization need only be computed once.
Given the sparse data structure for the Cholesky factor L, as supplied by the symbolic factor, the numeric factorization produces the entries in L so that
PAPT = LLT
Here P is the permutation matrix determined by the minimum degree ordering.
The numerical factorization can be carried out in one of two ways. By default, the standard factorization is performed based on a sparse compressed storage scheme. This is fully described in George and Liu (1981). Optionally, a multifrontal technique can be used. The multifrontal method requires more storage but will be faster in certain cases. The multifrontal factorization is based on the routines in Liu (1987). For a detailed description of this method, see Liu (1990), also Duff and Reid (1983, 1984), Ashcraft (1987), Ashcraft et al. (1987), and Liu (1986, 1989).
If an application requires that several linear systems be solved where the coefficient matrix is the same but the right-hand sides change, the options IMSL_RETURN_NUMERIC_FACTOR and IMSL_SUPPLY_NUMERIC_FACTOR can be used to precompute the Cholesky factor. Then the IMSL_SOLVE_ONLY option can be used to efficiently solve all subsequent systems.
Given the numeric factorization, the solution x is obtained by the following calculations:
Ly1 = Pb
Lty2 = y1
x = Pty2
The permutation information, P, is carried in the numeric factor structure.
Examples
Example 1
As an example consider the 5 × 5 coefficient matrix:
Let xT = (5, 4, 3, 2, 1) so that Ax = (55, 83, 103, 97, 82)T. The number of nonzeros in the lower triangle of A is nz = 10. The sparse coordinate form for the lower triangle is given by the following:
|
row |
0 |
1 |
2 |
2 |
3 |
3 |
4 |
4 |
4 |
4 |
|
col |
0 |
1 |
0 |
2 |
2 |
3 |
0 |
1 |
3 |
4 |
|
val |
10 |
20 |
1 |
30 |
4 |
40 |
2 |
3 |
5 |
50 |
Since this representation is not unique, an equivalent form would be as follows:
|
row |
3 |
4 |
4 |
4 |
0 |
1 |
2 |
2 |
3 |
4 |
|
col |
3 |
0 |
1 |
3 |
0 |
1 |
0 |
2 |
2 |
4 |
|
val |
40 |
2 |
3 |
5 |
10 |
20 |
1 |
30 |
4 |
50 |
#include <imsl.h>
#include <stdlib.h>
int main()
{
Imsl_f_sparse_elem a[] = {0, 0, 10.0,
1, 1, 20.0,
2, 0, 1.0,
2, 2, 30.0,
3, 2, 4.0,
3, 3, 40.0,
4, 0, 2.0,
4, 1, 3.0,
4, 3, 5.0,
4, 4, 50.0};
float b[] = {55.0, 83.0, 103.0, 97.0, 82.0};
int n = 5;
int nz = 10;
float *x;
x = imsl_f_lin_sol_posdef_coordinate (n, nz, a, b, 0);
imsl_f_write_matrix ("solution", 1, n, x, 0);
imsl_free (x);
}
Output
solution
1 2 3 4 5
5 4 3 2 1
Example 2
In this example, set A = E(2500, 50). Then solve the system Ax = bl and return the numeric factorization resulting from that call. Then solve the system Ax = b2 using the numeric factorization just computed. The ratio of execution time is printed. Be aware that timing results are highly machine dependent.
#include <imsl.h>
int main()
{
Imsl_f_sparse_elem *a;
Imsl_f_numeric_factor numeric_factor;
float *b_1;
float *b_2;
float *x_1;
float *x_2;
int n;
int ic;
int nz;
double time_1;
double time_2;
ic = 50;
n = ic*ic;
/* Generate two right hand sides */
b_1 = imsl_f_random_uniform (n*sizeof(*b_1), 0);
b_2 = imsl_f_random_uniform (n*sizeof(*b_2), 0);
/* Build coefficient matrix a */
a = imsl_f_generate_test_coordinate (n, ic, &nz,
IMSL_SYMMETRIC_STORAGE,
0);
/* Now solve Ax_1 = b_1 and return the numeric
factorization */
time_1 = imsl_ctime ();
x_1 = imsl_f_lin_sol_posdef_coordinate (n, nz, a, b_1,
IMSL_RETURN_NUMERIC_FACTOR, &numeric_factor,
0);
time_1 = imsl_ctime () - time_1;
/* Now solve Ax_2 = b_2 given the numeric
factorization */
time_2 = imsl_ctime ();
x_2 = imsl_f_lin_sol_posdef_coordinate (n, nz, a, b_2,
IMSL_SUPPLY_NUMERIC_FACTOR, &numeric_factor,
IMSL_SOLVE_ONLY,
0);
time_2 = imsl_ctime () - time_2;
printf("time_2/time_1 = %lf\n", time_2/time_1);
}
Output
time_2/time_1 = 0.037037
