sparse_cholesky_smp (complex)

 

more...

more...

Computes the Cholesky factorization of a sparse Hermitian positive definite matrix A by an OpenMP parallelized supernodal algorithm and solves the sparse Hermitian positive definite system of linear equations Ax = b.

Synopsis

#include <imsl.h>

f_complex *imsl_c_sparse_cholesky_smp (int n, int nz, Imsl_c_sparse_elem a[], f_complex b[], …, 0)

void imsl_free_snodal_symbolic_factor (Imsl_snodal_symbolic_factor  *sym_factor)

void imsl_c_free_numeric_factor (Imsl_c_numeric_factor *num_factor)

The type d_complex functions are imsl_z_sparse_cholesky_smp and imsl_z_free_numeric_factor.

Required Arguments

int n (Input)
The order of the input matrix.

int nz (Input)
Number of nonzeros in the lower triangle of the matrix.

Imsl_c_sparse_elem a[] (Input)
An array of length nz containing the location and value of each nonzero entry in the lower triangle of the matrix.

f_complex b[] (Input)
An array of length n containing the right-hand side.

Return Value

A pointer to the solution x of the sparse Hermitian 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>

f_complex  *imsl_c_sparse_cholesky_smp (int n, int nz, Imsl_c_sparse_elem a[], f_complex b[],

IMSL_CSC_FORMAT, int col_ptr[], int row_ind[], f_complex values[],

IMSL_PREORDERING, int preorder,

IMSL_RETURN_SYMBOLIC_FACTOR, Imsl_snodal_symbolic_factor *sym_factor,

IMSL_SUPPLY_SYMBOLIC_FACTOR, Imsl_snodal_symbolic_factor *sym_factor,

IMSL_SYMBOLIC_FACTOR_ONLY,

IMSL_RETURN_NUMERIC_FACTOR, Imsl_c_numeric_factor *num_factor,

IMSL_SUPPLY_NUMERIC_FACTOR, Imsl_c_numeric_factor *num_factor,

IMSL_NUMERIC_FACTOR_ONLY,

IMSL_SOLVE_ONLY,

IMSL_SMALLEST_DIAGONAL_ELEMENT, float *smallest_element,

IMSL_LARGEST_DIAGONAL_ELEMENT, float *largest_element,

IMSL_NUM_NONZEROS_IN_FACTOR, int *num_nonzeros,

IMSL_RETURN_USER, f_complex x[],

0)

Optional Arguments

IMSL_CSC_FORMAT, int col_ptr[], int row_ind[], f_complex values[] (Input)
Accept the coefficient matrix in compressed sparse column (CSC) format, as describedin the Compressed Sparse Column (CSC) Format section of the “Introduction” chapter of this manual.

IMSL_PREORDERING, int preorder (Input)
The variant of the Minimum Degree Ordering (MDO) algorithm used in the preordering of matrix A:

preorder	Method
0	George and Liu’s Quotient Minimum Degree algorithm.
1	A variant of George and Liu’s Quotient Minimum Degree algorithm using a preprocessing phase and external degrees.

Default: preorder = 0.

IMSL_RETURN_SYMBOLIC_FACTOR, Imsl_snodal_symbolic_factor *sym_factor (Output)
A pointer to a structure of type Imsl_snodal_symbolic_factor containing, on return, the supernodal symbolic factorization of the input matrix. A detailed description of the Imsl_snodal_symbolic_factor structure is given in the following table:

Table 1 – Strucuture Imsl_snodal_symbolic_factor

Parameter	Data Type	Description
nzsub	int **	A pointer to an array containing the compressed row subscripts of the non-zero off-diagonal elements of the Cholesky factor.
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]].
maxsub	int	The number of elements in array *nzsub that are used as subscripts. Note that the size of *nzsub can be larger than maxsub.
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 a description of alnz, see the description section of optional argument IMSL_RETURN_NUMERIC_FACTOR). For column j of the factor matrix, the starting and stopping indices of *alnz are stored in (*xlnz)[j] and (*xlnz)[j+1] respectively.
maxlnz	int	The number of non-zero off-diagonal elements in the Cholesky factor.
perm	int **	A pointer to an array of length n containing the permutation vector.
invp	int **	A pointer to an array of length n containing the inverse permutation vector.
multifrontal_space	int	This variable is not used in the current implementation.
nsuper	int	The number of supernodes in the Cholesky factor.
snode	int **	A pointer to an array of length n. Element (*snode)[j] contains the number of the fundamental supernode to which column j belongs.
snode_ptr	int **	A pointer to an array of length nsuper+1 containing the start column of each supernode.
nleaves	int	The number of leaves in the postordered elimination tree of the symmetrically permuted input matrix A.
etree_leaves	int **	A pointer to an array of length nleaves+1 containing the leaves of the elimination tree.

To free the memory allocated within this structure, use function imsl_free_snodal_symbolic_factor.

IMSL_SUPPLY_SYMBOLIC_FACTOR, Imsl_snodal_symbolic_factor *sym_factor (Input)
A pointer to a structure of type Imsl_snodal_symbolic_factor. This structure contains the symbolic factorization of the input matrix computed by imsl_c_sparse_cholesky_smp with the IMSL_RETURN_SYMBOLIC_FACTOR option. The structure is described in the IMSL_RETURN_SYMBOLIC_FACTOR optional argument description.
To free the memory allocated within this structure, use function imsl_free_snodal_symbolic_factor.

IMSL_SYMBOLIC_FACTOR_ONLY, (Input)
Compute the symbolic factorization of the input matrix and return. The argument b is ignored.

IMSL_RETURN_NUMERIC_FACTOR, Imsl_c_numeric_factor *num_factor (Output)
A pointer to a structure of type Imsl_c_numeric_factor containing, on return, the numeric factorization of the input matrix. A detailed description of the Imsl_c_numeric_factor structure is given in the IMSL_RETURN_NUMERIC_FACTOR optional argument description of function imsl_c_lin_sol_posdef_coordinate (complex). To free the memory allocated within this structure, use function imsl_c_free_numeric_factor.

IMSL_SUPPLY_NUMERIC_FACTOR, Imsl_c_numeric_factor *num_factor (Input)
A pointer to a structure of type Imsl_c_numeric_factor. This structure contains the numeric factorization of the input matrix computed by imsl_c_sparse_cholesky_smp with the IMSL_RETURN_NUMERIC_FACTOR option. The structure is described in the IMSL_RETURN_NUMERIC_FACTOR optional argument description of function imsl_lin_sol_posdef_coordinate (complex).

To free the memory allocated within this structure, use function imsl_c_free_numeric_factor.

IMSL_NUMERIC_FACTOR_ONLY, (Input)
Compute the numeric factorization of the input matrix and return. The argument b is ignored.

IMSL_SOLVE_ONLY, (Input)
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_SMALLEST_DIAGONAL_ELEMENT, float *smallest_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_c_sparse_cholesky_smp.

IMSL_LARGEST_DIAGONAL_ELEMENT, float *largest_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_c_sparse_cholesky_smp.

IMSL_NUM_NONZEROS_IN_FACTOR, int *num_nonzeros (Output)
A pointer to a scalar containing the total number of nonzeros in the factor.

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

Description

The function imsl_c_sparse_cholesky_smp solves a system of linear algebraic equations having a sparse Hermitian 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 exploiting OpenMP parallelism 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 supernodal data structure for the Cholesky factor, L. This step only requires the “pattern” of the sparse coefficient matrix, i.e., the locations of the nonzero elements but not any of the elements themselves. Thus, the val field in the Imsl_c_sparse_elem structure is ignored. If an application generates different sparse Hermitian 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 needs 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 = LLH

Here P is the permutation matrix determined by the minimum degree ordering.

The numerical factorization is an implementation of a parallel supernodal algorithm that combines a left-looking and a right-looking column computation scheme. This algorithm is described in detail in Rauber et al. (1999).

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

LHy2 = y1

x = PTy2

The permutation information, P, is carried in the numeric factor structure Imsl_c_numeric_factor.

Examples

Example 1

As a simple example of default use, consider the following Hermitian positive definite matrix

 

 

We construct the solution xT = (1 + i, 2 + 2i, 3 + 3i) to the system Ax = b by setting b:=Ax = (‑2 + 2i, 5 +15i, 36 + 28i)T. The number of nonzeros in the lower triangle is nz = 5. The solution is computed and printed.

 

#include <imsl.h>

 

int main()

{

int n = 3, nz = 5;

f_complex b[] = {{-2.0, 2.0}, {5.0, 15.0}, {36.0, 28.0}};

f_complex *x = NULL;

Imsl_c_sparse_elem a[] = {0, 0, {2.0, 0.0},

1, 1, {4.0, 0.0},

2, 2, {10.0, 0.0},

1, 0, {-1.0, -1.0},

2, 1, {1.0, -2.0}};

 

x = imsl_c_sparse_cholesky_smp (n, nz, a, b, 0);

 

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

 

imsl_free (x);

}

Output

 

Solution, x, of Ax = b

1 ( 1, 1)

2 ( 2, 2)

3 ( 3, 3)

Example 2

This example shows how a symbolic factor can be re-used. Consider matrix A, a Hermitian positive definite matrix with value 6 on the diagonal and value -1- i on its lower codiagonal and the lower band at distance 50 from the diagonal. At first, the system Ax = b is solved and the symbolic factorization of A is returned. Then, the system Cy = d with C = A+4*I, I the identity matrix, is solved using the symbolic factorization just computed. This is possible because A and C have the same nonzero structure and therefore also the same symbolic factorization. The solution errors are printed.

 

#include <imsl.h>

#include <stdlib.h>

#include <stdio.h>

 

int main()

{

int n, ic, nz, i, index;

float error_1, error_2;

f_complex *b = NULL, *d = NULL, *x = NULL, *y = NULL;

f_complex *mod_vector = NULL;

Imsl_c_sparse_elem *a = NULL, *c = NULL;

Imsl_snodal_symbolic_factor symbolic_factor;

 

ic = 50;

n = ic * ic;

 

mod_vector = (f_complex*) malloc (n * sizeof(f_complex));

 

/* Build coefficient matrix A */

a = imsl_c_generate_test_coordinate (n, ic,

&nz,

IMSL_SYMMETRIC_STORAGE,

0);

 

/* Build coefficient matrix C */

c = (Imsl_c_sparse_elem *) malloc (nz * sizeof (Imsl_c_sparse_elem));

 

for (i=0; i<nz; i++)

c[i] = a[i];

 

for (i=0; i<n; i++)

{

c[i].val.re = 10.0;

c[i].val.im = 0.0;

}

 

/* Form right hand side b */

for (i = 0; i < n; i++)

{

mod_vector[i].re = (float) (i % 5);

mod_vector[i].im = 0.0;

}

b = (f_complex *) imsl_c_mat_mul_rect_coordinate ("A*x",

IMSL_A_MATRIX, n, n, nz, a,

IMSL_X_VECTOR, n, mod_vector,

IMSL_SYMMETRIC_STORAGE,

0);

 

/* Form right hand side d */

d = (f_complex *) imsl_c_mat_mul_rect_coordinate ("A*x",

IMSL_A_MATRIX, n, n, nz, c,

IMSL_X_VECTOR, n, mod_vector,

IMSL_SYMMETRIC_STORAGE,

0);

 

/* Solve Ax = b and return the symbolic factorization */

x = imsl_c_sparse_cholesky_smp (n, nz, a, b,

IMSL_RETURN_SYMBOLIC_FACTOR, &symbolic_factor,

0);

 

/* Compute error |x-mod_vector| */

error_1 = imsl_c_vector_norm (n, x,

IMSL_SECOND_VECTOR, mod_vector,

IMSL_INF_NORM, &index,

0);

 

/* Solve Cy = d given the symbolic factorization */

y = imsl_c_sparse_cholesky_smp (n, nz, c, d,

IMSL_SUPPLY_SYMBOLIC_FACTOR, &symbolic_factor,

0);

 

/* Compute error |y-mod_vector| */

error_2 = imsl_c_vector_norm (n, y,

IMSL_SECOND_VECTOR, mod_vector,

IMSL_INF_NORM, &index,

0);

 

printf ("Solution error |x - mod_vector| = %e\n", error_1);

printf ("Solution error |y - mod_vector| = %e\n", error_2);

 

/* Free allocated memory */

if (mod_vector) free(mod_vector);

if (a) imsl_free (a);

if (c) free (c);

if (b) imsl_free (b);

if (d) imsl_free (d);

if (y) imsl_free (y);

if (x) imsl_free (x);

imsl_free_snodal_symbolic_factor(&symbolic_factor);

}

 

Output

 

Solution error |x - mod_vector| = 2.885221e-006

Solution error |y - mod_vector| = 2.386146e-006

Example 3

In this example, set A = E (2500, 50). First solve the system Ax = b1 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 times is printed. Be aware that timing results are highly machine dependent.

 

#include <imsl.h>

#include <stdio.h>

#include <omp.h>

 

int main()

{

int n, ic, nz, i, index;

float *rand_vec = NULL;

double time_1, time_2;

f_complex b_1[2500], b_2[2500], *x_1 = NULL, *x_2 = NULL;

Imsl_c_sparse_elem *a = NULL;

Imsl_c_numeric_factor numeric_factor;

 

ic = 50;

n = ic * ic;

index = 0;

 

/* Generate two right hand sides */

imsl_random_seed_set (1234567);

rand_vec = imsl_f_random_uniform (4 * n, 0);

 

for (i = 0; i < n; i++) {

b_1[i].re = rand_vec[index++];

b_1[i].im = rand_vec[index++];

b_2[i].re = rand_vec[index++];

b_2[i].im = rand_vec[index++];

}

 

/* Build coefficient matrix a */

a = imsl_c_generate_test_coordinate (n, ic, &nz,

IMSL_SYMMETRIC_STORAGE,

0);

 

/* Now solve Ax_1 = b_1 and return the numeric factorization */

time_1 = omp_get_wtime();

 

x_1 = imsl_c_sparse_cholesky_smp (n, nz, a, b_1,

IMSL_RETURN_NUMERIC_FACTOR, &numeric_factor,

0);

 

time_1 = omp_get_wtime() - time_1;

 

/* Now solve Ax_2 = b_2 given the numeric factorization */

time_2 = omp_get_wtime();

 

x_2 = imsl_c_sparse_cholesky_smp (n, nz, a, b_2,

IMSL_SUPPLY_NUMERIC_FACTOR, &numeric_factor,

IMSL_SOLVE_ONLY,

0);

 

time_2 = omp_get_wtime() - time_2;

 

printf("time_2/time_1 = %lf\n", time_2/time_1);

 

/* Free memory */

if (rand_vec) imsl_free(rand_vec);

if (x_1) imsl_free(x_1);

if (x_2) imsl_free(x_2);

if (a) imsl_free(a);

imsl_c_free_numeric_factor(&numeric_factor);

}

 

Output

 

time_2/time_1 = 0.025771

 

Fatal Errors

IMSL_BAD_SQUARE_ROOT

A zero or negative square root has occurred during the factorization. The coefficient matrix is not positive definite.