Computes the LU factorization of a general sparse matrix by a left-looking column method using OpenMP parallelism, and solves the real sparse linear system of equations Ax = b .

The type double functions are imsl_d_superlu_smp and imsl_d_superlu_smp_factor_free.

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.

IMSL_TRANSPOSE, int transpose (Input)
Indicates if the transposed problem ATx = b is to be solved. This option can be used in conjunction with either of the options that supply the factorization.

IMSL_SNODE_PREDICTION, int snode_prediction (Input)
Indicates which scheme is used to predict the number of nonzeros in the L supernodes.

The steps imsl_f_superlu_smp uses to solve linear systems are identical to the steps described in the documentation of the serial version imsl_f_superlu.

Function imsl_f_superlu_smp uses a supernodal storage scheme for the LU factorization of matrix A. In contrast to the sequential version, the consecutive columns and supernodes of the L and U factors might not be stored contiguously in memory. Thus, in addition to the pointers to the beginning of each column or supernode, also pointers to the end of each column or supernode are needed. The factorization is contained in structure Imsl_f_super_lu_smp_factor and its two sub-structures Imsl_f_hbp_format and Imsl_f_scp_format. Following is a short description of these structures:

Structure Imsl_d_super_lu_smp_factor and its two sub-structures are defined similarly by replacing float with double, Imsl_f_hbp_format with Imsl_d_hbp_format, and Imsl_f_scp_format with Imsl_d_scp_format in their respective definitions.

In contrast to the sequential version, the numerical factorization phase of the LU decomposition is parallelized. Since a dynamic memory expansion as in the serial case is difficult to implement for the parallel code, the estimated sizes of array rowind for the L and of arrays rowind and nzval for the U factor (see structures Imsl_f_scp_format and Imsl_f_hbp_format above) must be predetermined by the user via elements 6 and 7 of the performance tuning array sp_ienv.

In order to ensure that the columns of each L supernode are stored contiguously in memory, a static or dynamic prediction scheme for the size of the L supernodes can be used. The static version, which function imsl_f_superlu_smp uses by default, exploits the observation that for any row permutation P in PA = LU, the nonzero structure of L is contained in that of the Householder matrix H from the Householder sparse QR factorization A = QR. Furthermore, it can be shown that each fundamental supernode in L is always contained in a fundamental supernode of H. Therefore, the storage requirement for the L supernodes and array nzval in the L factor respectively can be estimated and allocated prior to the factorization based on the size of the H supernodes. The algorithm used to compute the supernode partition and the size of the supernodes in H is almost linear in the number of nonzeros of matrix A.

In practice, the above static prediction scheme is quite tight for most problems. However, if the number of nonzeros in H greatly exceeds the number of nonzeros in L, the user can try a dynamic prediction scheme by setting optional argument IMSL_SNODE_PREDICTION to 1. This scheme still uses the supernode partition in H, but dynamically searches the supernodal graph of L to obtain a much tighter upper bound for the required storage. Use of the dynamic scheme requires the user to define the size of array nzval in the L factor via element 5 of the performance tuning array sp_ienv.

For a complete description of the parallel algorithm, see Demmel et al. (1999c).

Copyright (c) 2003, The Regents of the University of California, through Lawrence Berkeley National Laboratory (subject to receipt of any required approvals from U.S. Dept. of Energy).

All rights reserved.

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

(1) Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.

(2) Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.

(3) Neither the name of Lawrence Berkeley National Laboratory, U.S. Dept. of Energy nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

The LU factorization of the sparse 6×6 matrix

is computed.

Let y = (1,2,3,4,5,6)T, so that b1 := Ay = (10,7,45,33,-34,31)T and b2 := ATy = (-9,8,39,13,1,21)T.

The LU factorization of A is used to solve the sparse linear systems Ax = b1 and ATx = b2.

This example uses the matrix A = E(1000,10) to show how the LU factorization of A can be used to solve a linear system with the same coefficient matrix A but different right-hand side. Maximum absolute errors are printed. After the computations, the space allocated for the LU factorization is freed via function imsl_f_superlu_smp_factor_free.