Engineering and Scientific Subroutine Library for AIX Version 3 Release 3: Guide and Reference

DSKFS--Symmetric Sparse Matrix Factorization, Determinant, and Solve Using Skyline Storage Mode

This subroutine can perform either or both of the following functions for symmetric sparse matrix A, stored in skyline storage mode, and for vectors x and b:

• Factor A and, optionally, compute the determinant of A.
• Solve the system Ax = b using the results of the factorization of matrix A, produced on this call or a preceding call to this subroutine.

You have the choice of using either Gaussian elimination or Cholesky decomposition. You also have the choice of using profile-in or diagonal-out skyline storage mode for A on input or output.

Note:

The input to the solve performed by this subroutine must be the output from the factorization performed by this subroutine.

Syntax

On Entry

n

is the order of symmetric sparse matrix A. Specified as: a fullword integer; n >= 0.

a

is the array, referred to as A, containing one of three forms of the upper triangular part of symmetric sparse matrix A, depending on the type of computation performed, where:

• If you are doing a factor and solve or a factor only, and if IPARM(3) = 0, then A contains the unfactored upper triangle of symmetric sparse matrix A.
• If you are doing a factor only, and if IPARM(3) > 0, then A contains the partially factored upper triangle of symmetric sparse matrix A. The first IPARM(3) columns in the upper triangle of A are already factored. The remaining columns are factored in this computation.
• If you are doing a solve only, then A contains the factored upper triangle of sparse matrix A, produced by a preceding call to this subroutine.

In each case:

If IPARM(4) = 0, diagonal-out skyline storage mode is used for A.

If IPARM(4) = 1, profile-in skyline storage mode is used for A.

Specified as: a one-dimensional array of (at least) length na, containing long-precision real numbers.

na

is the length of array A. Specified as: a fullword integer; na >= 0 and na >= (IDIAG(n+1)-1).

idiag

is the array, referred to as IDIAG, containing the relative positions of the diagonal elements of matrix A (in one of its three forms) in array A. Specified as: a one-dimensional array of (at least) length n+1, containing fullword integers.

iparm

is an array of parameters, IPARM(i), where:

• IPARM(1) indicates whether certain default values for iparm and rparm are used by this subroutine, where:

  If IPARM(1) = 0, the following default values are used. For restrictions, see Notes.

  IPARM(2) = 0

  IPARM(3) = 0

  IPARM(4) = 0

  IPARM(5) = 0

  IPARM(10) = 0

  IPARM(11) = -1

  IPARM(12) = -1

  IPARM(13) = -1

  IPARM(14) = -1

  IPARM(15) = 0

  RPARM(10) = 10^-12

  If IPARM(1) = 1, the default values are not used.

• IPARM(2) indicates the type of computation performed by this subroutine. The following table gives the IPARM(2) values for each variation:
  

  Type of Computation	Gaussian Elimination Ax = b	Gaussian Elimination Ax = b and Determinant(A)	Cholesky Decomposition Ax = b	Cholesky Decomposition Ax = b and Determinant(A)
  Factor and Solve	0	10	100	110
  Factor Only	1	11	101	111
  Solve Only	2	N/A	102	N/A

• IPARM(3) indicates whether a full or partial factorization is performed on matrix A, where:

  If IPARM(3) = 0, and:

  If you are doing a factor and solve or a factor only, then a full factorization is performed for matrix A on rows and columns 1 through n.

  If you are doing a solve only, this argument has no effect on the computation, but must be set to 0.

  If IPARM(3) > 0, and you are doing a factor only, then a partial factorization is performed on matrix A. Rows 1 through IPARM(3) of columns 1 through IPARM(3) in matrix A must be in factored form from a preceding call to this subroutine. The factorization is performed on rows IPARM(3)+1 through n and columns IPARM(3)+1 through n. For an illustration, see Notes.

• IPARM(4) indicates the input storage mode used for matrix A. This determines the arrangement of data in arrays A and IDIAG on input, where:

  If IPARM(4) = 0, diagonal-out skyline storage mode is used.

  If IPARM(4) = 1, profile-in skyline storage mode is used.

• IPARM(5) indicates the output storage mode used for matrix A. This determines the arrangement of data in arrays A and IDAIG on output, where:

  If IPARM(5) = 0, diagonal-out skyline storage mode is used.

  If IPARM(5) = 1, profile-in skyline storage mode is used.

• IPARM(6) through IPARM(9) are reserved.
• IPARM(10) has the following meaning, where:

  If you are doing a factor and solve or a factor only, then IPARM(10) indicates whether certain default values for iparm and rparm are used by this subroutine, where:

  If IPARM(10) = 0, the following default values are used. For restrictions, see Notes.

  IPARM(11) = -1

  IPARM(12) = -1

  IPARM(13) = -1

  IPARM(14) = -1

  IPARM(15) = 0

  RPARM(10) = 10^-12

  If IPARM(10) = 1, the default values are not used.

  If you are doing a solve only, this argument is not used.

• IPARM(11) through IPARM(15) have the following meaning, where:

  If you are doing a factor and solve or a factor only, then IPARM(11) through IPARM(15) control the type of processing to apply to pivot elements occurring in regions 1 through 5, respectively. The pivot elements are d_kk for Gaussian elimination and r_kk for Cholesky decomposition for k = 1, n when doing a full factorization, and they are k = IPARM(3)+1, n when doing a partial factorization. The region in which a pivot element falls depends on the sign and magnitude of the pivot element. The regions are determined by RPARM(10). For a description of the regions and associated pivot values, see Notes. For each region i for i = 1,5, where the pivot occurs in region i, the processing applied to the pivot element is determined by IPARM(10+i), where:

  If IPARM(10+i) = -1, the pivot element is trapped and computational error 2126 is generated. See Error Conditions.

  If IPARM(10+i) = 0, processing continues normally.

  Note:

  A value of 0 is not permitted for region 3, because if processing continues, a divide-by-zero exception occurs. In addition, if you are doing a Cholesky decomposition, a value of 0 is not permitted in regions 1 and 2, because a square root exception occurs.

  If IPARM(10+i) = 1, the pivot element is replaced with the value in RPARM(10+i), and processing continues normally.

  If you are doing a solve only, these arguments are not used.

• IPARM(16) through IPARM(25), see On Return.

Specified as: a one-dimensional array of (at least) length 25, containing fullword integers, where:

IPARM(1) = 0 or 1

IPARM(2) = 0, 1, 2, 10, 11, 100, 101, 102, 110, or 111

If IPARM(2) = 0, 2, 10, 100, 102, or 110, then IPARM(3) = 0

If IPARM(2) = 1, 11, 101, or 111, then 0 <= IPARM(3) <= n

IPARM(4), IPARM(5) = 0 or 1

If IPARM(2) = 0, 1, 10, or 11, then:

IPARM(10) = 0 or 1

IPARM(11), IPARM(12) = -1, 0, or 1

IPARM(13) = -1 or 1

IPARM(14), IPARM(15) = -1, 0, or 1

If IPARM(2) = 100, 101, 110, or 111, then:

IPARM(10) = 0 or 1

IPARM(11), IPARM(12), IPARM(13) = -1 or 1

IPARM(14), IPARM(15) = -1, 0, or 1

rparm

is an array of parameters, RPARM(i), where:

• RPARM(1) through RPARM(9) are reserved.
• RPARM(10) has the following meaning, where:

  If you are doing a factor and solve or a factor only, RPARM(10) is the tolerance value for small pivots. This sets the bounds for the pivot regions, where pivots are processed according to the options you specify for the five regions in IPARM(11) through IPARM(15), respectively. The suggested value is 10^-15 <= IPARM(10) <= 1.

  If you are doing a solve only, this argument is not used.

• RPARM(11) through RPARM(15) have the following meaning, where:

  If you are doing a factor and solve or a factor only, RPARM(11) through RPARM(15) are the fix-up values to use for the pivots in regions 1 through 5, respectively. For each RPARM(10+i) for i = 1,5, where the pivot occurs in region i:

  If IPARM(10+i) = 1, the pivot is replaced with RPARM(10+i), where |RPARM(10+i)| should be a sufficiently large nonzero value to avoid overflow when calculating the reciprocal of the pivot. For Gaussian elimination, the suggested value is 10^-15 <= |RPARM(10+i)| <= 1. For Cholesky decomposition, the value must be RPARM(10+i) > 0.

  If IPARM(10+i) <> 1, RPARM(10+i) is not used.

  If you are doing a solve only, these arguments are not used.

• RPARM(16) through RPARM(25), see On Return.

Specified as: a one-dimensional array of (at least) length 25, containing long-precision real numbers, where if IPARM(2) = 0, 1, 10, 11, 100, 101, 110, or 111, then:

RPARM(10) >= 0.0

If IPARM(2) = 0, 1, 10, or 11, then RPARM(11) through RPARM(15) <> 0.0

If IPARM(2) = 100, 101, 110, or 111, then RPARM(11) through RPARM(15) > 0.0

aux

has the following meaning:

If naux = 0 and error 2015 is unrecoverable, aux is ignored.

Otherwise, it is the storage work area used by this subroutine. Its size is specified by naux.

Specified as: an area of storage, containing long-precision real numbers.

naux

is the size of the work area specified by aux--that is, the number of elements in aux. Specified as: a fullword integer, where:

If naux = 0 and error 2015 is unrecoverable, DSKFS dynamically allocates the work area used by this subroutine. The work area is deallocated before control is returned to the calling program.

Otherwise, If you are doing a factor only, you can use naux >= n; however, for optimal performance, use naux >= 3n.

If you are doing a factor and solve or a solve only, use naux >= 3n+4mbx.

For further details on error handling and the special factor-only case, see Notes.

bx

has the following meaning, where:

If you are doing a factor and solve or a solve only, bx is the array, containing the mbx right-hand side vectors b of the system Ax = b. Each vector b is length n and is stored in the corresponding column of the array.

If you are doing a factor only, this argument is not used in the computation.

Specified as: an ldbx by (at least) mbx array, containing long-precision real numbers.

ldbx

has the following meaning, where:

If you are doing a factor and solve or a solve only, ldbx is the leading dimension of the array specified for bx.

If you are doing a factor only, this argument is not used in the computation.

Specified as: a fullword integer; ldbx >= n and:

If mbx <> 0, then ldbx > 0.

If mbx = 0, then ldbx >= 0.

mbx

has the following meaning, where:

If you are doing a factor and solve or a solve only, mbx is the number of right-hand side vectors, b, in the array specified for bx.

If you are doing a factor only, this argument is not used in the computation.

Specified as: a fullword integer; mbx >= 0.

On Return

a

is the array, referred to as A, containing the upper triangular part of symmetric sparse matrix A in LDL^T or R^TR factored form, where:

If IPARM(5) = 0, diagonal-out skyline storage mode is used for A.

If IPARM(5) = 1, profile-in skyline storage mode is used for A.

(If mbx = 0 and you are doing a solve only, then a is unchanged on output.) Returned as: a one-dimensional array of (at least) length na, containing long-precision real numbers.

idiag

is the array, referred to as IDIAG, containing the relative positions of the diagonal elements of the factored output matrix A in array A. (If mbx = 0 and you are doing a solve only, then idiag is unchanged on output.)

Returned as: a one-dimensional array of (at least) length n+1, containing fullword integers.

iparm

is an array of parameters, IPARM(i), where:

• IPARM(1) through IPARM(15) are unchanged.
• IPARM(16) has the following meaning, where:

  If you are doing a factor and solve or a factor only, and:

  If IPARM(16) = -1, your factorization did not complete successfully, resulting in computational error 2126.

  If IPARM(16) > 0, it is the row number k, in which the maximum absolute value of the ratio a_kk/d_kk for Gaussian elimination and a_kk/r_kk for Cholesky decomposition occurred, where:

  If IPARM(3) = 0, k can be any of the rows, 1 through n, in the full factorization.

  If IPARM(3) > 0, k can be any of the rows, IPARM(3)+1 through n, in the partial factorization.

  If you are doing a solve only, this argument is not used in the computation and is unchanged.

• IPARM(17) through IPARM(20) are reserved.
• IPARM(21) through IPARM(25) have the following meaning, where:

  If you are doing a factor and solve or a factor only, IPARM(21) through IPARM(25) have the following meanings for each region i for i = 1,5, respectively:

  If IPARM(20+i) = -1, your factorization did not complete successfully, resulting in computational error 2126.

  If IPARM(20+i) >= 0, it is the number of pivots in region i for the columns that were factored in matrix A, where:

  If IPARM(3) = 0, columns 1 through n were factored in the full factorization.

  If IPARM(3) > 0, columns IPARM(3)+1 through n were factored in the partial factorization.

  If you are doing a solve only, these arguments are not used in the computation and are unchanged.

Returned as: a one-dimensional array of (at least) length 25, containing fullword integers.

rparm

is an array of parameters, RPARM(i), where:

• RPARM(1) through RPARM(15) are unchanged.
• RPARM(16) has the following meaning, where:

  If you are doing a factor and solve or a factor only, and:

  If RPARM(16) = 0.0, your factorization did not complete successfully, resulting in computational error 2126.

  If |RPARM(16)| > 0.0, it is the ratio for row k, a_kk/d_kk for Gaussian elimination and a_kk/r_kk for Cholesky decomposition, having the maximum absolute value. Row k is indicated in IPARM(16), and:

  If IPARM(3) = 0, the ratio corresponds to one of the rows, 1 through n, in the full factorization.

  If IPARM(3) > 0, the ratio corresponds to one of the rows, IPARM(3)+1 through n, in the partial factorization.

  If you are doing a solve only, this argument is not used in the computation and is unchanged.

• RPARM(17) and RPARM(18) have the following meaning, where:

  If you are computing the determinant of matrix A, then RPARM(17) is the mantissa, detbas, and RPARM(18) is the power of 10, detpwr, used to express the value of the determinant: detbas(10^detpwr), where 1 <= detbas < 10. Also:

  If IPARM(3) = 0, the determinant is computed for columns 1 through n in the full factorization.

  If IPARM(3) > 0, the determinant is computed for columns IPARM(3)+1 through n in the partial factorization.

  If you are not computing the determinant of matrix A, these arguments are not used in the computation and are unchanged.

• RPARM(19) through RPARM(25) are reserved.

Returned as: a one-dimensional array of (at least) length 25, containing long-precision real numbers.

bx

has the following meaning, where:

If you are doing a factor and solve or a solve only, bx is the array, containing the mbx solution vectors x of the system Ax = b. Each vector x is length n and is stored in the corresponding column of the array. (If mbx = 0, then bx is unchanged on output.)

If you are doing a factor only, this argument is not used in the computation and is unchanged.

Returned as: an ldbx by (at least) mbx array, containing long-precision real numbers.

Notes

1. When doing a solve only, you should specify the same factorization method in IPARM(2), Gaussian elimination or Cholesky decomposition, that you specified for your factorization on a previous call to this subroutine.
2. If you set either IPARM(1) = 0 or IPARM(10) = 0, indicating you want to use the default values for IPARM(11) through IPARM(15) and RPARM(10), then:
   • Matrix A must be positive definite.
   • No pivots are fixed, using RPARM(11) through RPARM(15) values.
   • No small pivots are tolerated; that is, the value should be |pivot| > RPARM(10).
3. Many of the input and output parameters for iparm and rparm are defined for the five pivot regions handled by this subroutine. The limits of the regions are based on RPARM(10), as shown in Figure 12. The pivot values in each region are:

   Region 1: pivot < -RPARM(10)

   Region 2: -RPARM(10) <= pivot < 0

   Region 3: pivot = 0

   Region 4: 0 < pivot <= RPARM(10)

   Region 5: pivot > RPARM(10)

   Figure 12. Five Pivot Regions
   

4. The IPARM(4) and IPARM(5) arguments allow you to specify the same or different skyline storage modes for your input and output arrays for matrix A. This allows you to change storage modes as needed. However, if you are concerned with performance, you should use diagonal-out skyline storage mode for both input and output, if possible, because there is less overhead.

   For a description of how sparse matrices are stored in skyline storage mode, see Profile-In Skyline Storage Mode and Diagonal-Out Skyline Storage Mode. Those descriptions use different array and variable names from the ones used here. To relate the two sets, use the following table:
   

   Name Here	Name in the Storage Description
   A	AU
   na	nu
   IDIAG	IDU

5. Following is an illustration of the portion of matrix A factored in the partial factorization when IPARM(3) > 0. In this case, the subroutine assumes that rows and columns 1 through IPARM(3) are already factored and that rows and columns IPARM(3)+1 through n are to be factored in this computation.

   
   
   

   
   

   You use the partial factorization function when, for design or storage reasons, you must factor the matrix A in stages. When doing a partial factorization, you must use the same skyline storage mode for all parts of the matrix as it is progressively factored.

6. Your various arrays must have no common elements; otherwise, results are unpredictable.
7. You have the option of having the minimum required value for naux dynamically returned to your program. For details, see Using Auxiliary Storage in ESSL.

Function

This subroutine can factor, compute the determinant of, and solve symmetric sparse matrix A, stored in skyline storage mode. It can use either Gaussian elimination or Cholesky decomposition. For all computations, input matrix A can be stored in either diagonal-out or profile-in skyline storage mode. Output matrix A can also be stored in either of these modes and can be different from the mode used for input.

For Gaussian elimination, matrix A is factored into the following form using specified pivot processing:

A = LDL^T

where:

D is a diagonal matrix.

L is a lower triangular matrix.

The transformed matrix A, factored into its LDL^T form, is stored in packed format in array A, such that the inverse of the diagonal matrix D is stored in the corresponding elements of array A. The off-diagonal elements of the unit upper triangular matrix L^T are stored in the corresponding off-diagonal elements of array A.

For Cholesky decomposition, matrix A is factored into the following form using specified pivot processing:

A = R^TR

where R is an upper triangular matrix

The transformed matrix A, factored into its R^TR form, is stored in packed format in array A, such that the inverse of the diagonal elements of the upper triangular matrix R is stored in the corresponding elements of array A. The off-diagonal elements of matrix R are stored in the corresponding off-diagonal elements of array A.

The partial factorization of matrix A, which you can do when you specify the factor-only option, assumes that the first IPARM(3) rows and columns are already factored in the input matrix. It factors the remaining n-IPARM(3) rows and columns in matrix A. (See Notes for an illustration.) It updates only the elements in array A corresponding to the part of matrix A that is factored.

The determinant can be computed with any of the factorization computations. With a full factorization, you get the determinant for the whole matrix. With a partial factorization, you get the determinant for only that part of the matrix factored in this computation.

The system Ax = b, having multiple right-hand sides, is solved for x using the transformed matrix A produced by this call or a subsequent call to this subroutine.

See references [9], [12], [25], [47], [71]. If n is 0, no computation is performed. If mbx is 0, no solve is performed.

Error Conditions

Resource Errors

• Error 2015 is unrecoverable, naux = 0, and unable to allocate work area.
• Unable to allocate internal work area.

Computational Errors

1. If a pivot occurs in region i for i = 1,5 and IPARM(10+i) = 1, the pivot value is replaced with RPARM(10+i), an attention message is issued, and processing continues.
2. Unacceptable pivot values occurred in the factorization of matrix A.
   • One or more diagonal elements of D or R contains unacceptable pivots and no valid fixup is applicable. The row number i of the first unacceptable pivot element is identified in the computational error message.
   • The return code is set to 2.
   • i can be determined at run time by use of the ESSL error-handling facilities. To obtain this information, you must use ERRSET to change the number of allowable errors for error code 2126 in the ESSL error option table; otherwise, the default value causes your program to terminate when this error occurs. For details, see What Can You Do about ESSL Computational Errors?.

Input-Argument Errors

1. n < 0
2. na < 0
3. IDIAG(n+1) > na+1
4. IDIAG(i+1) <= IDIAG(i) for i = 1, n
5. IDIAG(i+1) > IDIAG(i)+i and IPARM(4) = 0 for i = 1, n
6. IDIAG(i) > IDIAG(i-1)+i and IPARM(4) = 1 for i = 2, n
7. IPARM(1) <> 0 or 1
8. IPARM(2) <> 0, 1, 2, 10, 11, 100, 101, 102, 110, or 111
9. IPARM(3) < 0
10. IPARM(3) > n
11. IPARM(3) > 0 and IPARM(2) <> 1, 11, 101, or 111
12. IPARM(4), IPARM(5) <> 0 or 1
13. IPARM(2) = 0, 1, 10, or 11 and:

    IPARM(10) <> 0 or 1

    IPARM(11), IPARM(12) <> -1, 0, or 1

    IPARM(13) <> -1 or 1

    IPARM(14), IPARM(15) <> -1, 0, or 1

    RPARM(10) < 0.0

    RPARM(10+i) = 0.0 and IPARM(10+i) = 1 for i = 1,5

14. IPARM(2) = 100, 101, 110, or 111 and:

    IPARM(10) <> 0 or 1

    IPARM(11), IPARM(12), IPARM(13) <> -1 or 1

    IPARM(14), IPARM(15) <> -1, 0, or 1

    RPARM(10) < 0.0

    RPARM(10+i) <= 0.0 and IPARM(10+i) = 1 for i = 1,5

15. IPARM(2) = 0, 2, 10, 100, 102, or 110 and:

    ldbx <= 0 and mbx <> 0 and n <> 0

    ldbx < 0 and mbx = 0

    ldbx < n and mbx <> 0

    mbx < 0

16. Error 2015 is recoverable or naux<>0, and naux is too small--that is, less than the minimum required value. Return code 1 is returned if error 2015 is recoverable.

Example 1

This example shows how to factor a 9 by 9 symmetric sparse matrix A and solve the system Ax = b with three right-hand sides. It uses Gaussian elimination. The default values are used for IPARM and RPARM. Input matrix A, shown here, is stored in diagonal-out skyline storage mode. Matrix A is:

            *                                             *
            | 1.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0  0.0 |
            | 1.0  2.0  2.0  2.0  1.0  1.0  0.0  1.0  0.0 |
            | 1.0  2.0  3.0  3.0  2.0  2.0  0.0  2.0  0.0 |
            | 1.0  2.0  3.0  4.0  3.0  3.0  0.0  3.0  0.0 |
            | 0.0  1.0  2.0  3.0  4.0  4.0  1.0  4.0  0.0 |
            | 0.0  1.0  2.0  3.0  4.0  5.0  2.0  5.0  1.0 |
            | 0.0  0.0  0.0  0.0  1.0  2.0  3.0  3.0  2.0 |
            | 0.0  1.0  2.0  3.0  4.0  5.0  3.0  7.0  3.0 |
            | 0.0  0.0  0.0  0.0  0.0  1.0  2.0  3.0  4.0 |
            *                                             *

Output matrix A, shown here, is in LDL^T factored form with D^-1 on the diagonal, and is stored in diagonal-out skyline storage mode. Matrix A is:

             *                                             *
             | 1.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0  0.0 |
             | 1.0  1.0  1.0  1.0  1.0  1.0  0.0  1.0  0.0 |
             | 1.0  1.0  1.0  1.0  1.0  1.0  0.0  1.0  0.0 |
             | 1.0  1.0  1.0  1.0  1.0  1.0  0.0  1.0  0.0 |
             | 0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  0.0 |
             | 0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0 |
             | 0.0  0.0  0.0  0.0  1.0  1.0  1.0  1.0  1.0 |
             | 0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0 |
             | 0.0  0.0  0.0  0.0  0.0  1.0  1.0  1.0  1.0 |
             *                                             *

Call Statement and Input

            N  A  NA  IDIAG  IPARM  RPARM  AUX NAUX  BX  LDBX MBX
            |  |  |     |      |      |     |   |    |    |    |
CALL DSKFS( 9, A, 33, IDIAG, IPARM, RPARM, AUX, 39 , BX , 12 , 3 )

A        =  (1.0, 2.0, 1.0, 3.0, 2.0, 1.0, 4.0, 3.0, 2.0, 1.0, 4.0,
             3.0, 2.0, 1.0, 5.0, 4.0, 3.0, 2.0, 1.0, 3.0, 2.0, 1.0,
             7.0, 3.0, 5.0, 4.0, 3.0, 2.0, 1.0, 4.0, 3.0, 2.0, 1.0)
IDIAG    =  (1, 2, 4, 7, 11, 15, 20, 23, 30, 34)
IPARM    =  (0, . , . , . , . , . , . , . , . , . , . , . , . , . ,
             . , . , . , . , . , . , . , . , . , . , . )

        *                     *
        |  4.00   8.00  12.00 |
        | 10.00  20.00  30.00 |
        | 15.00  30.00  45.00 |
        | 19.00  38.00  57.00 |
        | 19.00  38.00  57.00 |
BX   =  | 23.00  46.00  69.00 |
        | 11.00  22.00  33.00 |
        | 28.00  56.00  84.00 |
        | 10.00  20.00  30.00 |
        |   .      .      .   |
        |   .      .      .   |
        |   .      .      .   |
        *                     *

Output

        *                  *
        | 1.00  2.00  3.00 |
        | 1.00  2.00  3.00 |
        | 1.00  2.00  3.00 |
        | 1.00  2.00  3.00 |
        | 1.00  2.00  3.00 |
BX   =  | 1.00  2.00  3.00 |
        | 1.00  2.00  3.00 |
        | 1.00  2.00  3.00 |
        | 1.00  2.00  3.00 |
        |  .     .     .   |
        |  .     .     .   |
        |  .     .     .   |
        *                  *

Example 2

This example shows how to factor the 9 by 9 symmetric sparse matrix A from Example 1, solve the system Ax = b with three right-hand sides, and compute the determinant of A. It uses Gaussian elimination. The default values for pivot processing are used for IPARM. Input matrix A is stored in profile-in skyline storage mode. Output matrix A is in LDL^T factored form with D^-1 on the diagonal, and is stored in diagonal-out skyline storage mode. It is the same as output matrix A in Example 1.

Call Statement and Input

            N  A  NA  IDIAG  IPARM  RPARM  AUX NAUX  BX  LDBX MBX
            |  |  |     |      |      |     |   |    |    |    |
CALL DSKFS( 9, A, 33, IDIAG, IPARM, RPARM, AUX, 39 , BX , 12 , 3 )

A        =  (1.0, 1.0, 2.0, 1.0, 2.0, 3.0, 1.0, 2.0, 3.0, 4.0, 1.0,
             2.0, 3.0, 4.0, 1.0, 2.0, 3.0, 4.0, 5.0, 1.0, 2.0, 3.0,
             1.0, 2.0, 3.0, 4.0, 5.0, 3.0, 7.0, 1.0, 2.0, 3.0, 4.0)
IDIAG    =  (1, 3, 6, 10, 14, 19, 22, 29, 33, 34)
IPARM    =  (1, 10, 0, 1, 0, . , . , . , . , 0, . , . , . , . , . ,
             . , . , . , . , . , . , . , . , . , . )

        *                     *
        |  4.00   8.00  12.00 |
        | 10.00  20.00  30.00 |
        | 15.00  30.00  45.00 |
        | 19.00  38.00  57.00 |
        | 19.00  38.00  57.00 |
BX   =  | 23.00  46.00  69.00 |
        | 11.00  22.00  33.00 |
        | 28.00  56.00  84.00 |
        | 10.00  20.00  30.00 |
        |   .      .      .   |
        |   .      .      .   |
        |   .      .      .   |
        *                     *

Output

Example 3

This example shows how to factor a 9 by 9 negative-definite symmetric sparse matrix A, solve the system Ax = b with three right-hand sides, and compute the determinant of A. It uses Gaussian elimination. (Default values for pivot processing are not used for IPARM because A is negative-definite.) Input matrix A, shown here, is stored in diagonal-out skyline storage mode. Matrix A is:

         *                                                      *
         | -1.0  -1.0  -1.0  -1.0   0.0   0.0   0.0   0.0   0.0 |
         | -1.0  -2.0  -2.0  -2.0  -1.0  -1.0   0.0  -1.0   0.0 |
         | -1.0  -2.0  -3.0  -3.0  -2.0  -2.0   0.0  -2.0   0.0 |
         | -1.0  -2.0  -3.0  -4.0  -3.0  -3.0   0.0  -3.0   0.0 |
         |  0.0  -1.0  -2.0  -3.0  -4.0  -4.0  -1.0  -4.0   0.0 |
         |  0.0  -1.0  -2.0  -3.0  -4.0  -5.0  -2.0  -5.0  -1.0 |
         |  0.0   0.0   0.0   0.0  -1.0  -2.0  -3.0  -3.0  -2.0 |
         |  0.0  -1.0  -2.0  -3.0  -4.0  -5.0  -3.0  -7.0  -3.0 |
         |  0.0   0.0   0.0   0.0   0.0  -1.0  -2.0  -3.0  -4.0 |
         *                                                      *

Output matrix A, shown here, is in LDL^T factored form with D^-1 on the diagonal, and is stored in diagonal-out skyline storage mode. Matrix A is:

         *                                                      *
         | -1.0   1.0   1.0   1.0   0.0   0.0   0.0   0.0   0.0 |
         |  1.0  -1.0   1.0   1.0   1.0   1.0   0.0   1.0   0.0 |
         |  1.0   1.0  -1.0   1.0   1.0   1.0   0.0   1.0   0.0 |
         |  1.0   1.0   1.0  -1.0   1.0   1.0   0.0   1.0   0.0 |
         |  0.0   1.0   1.0   1.0  -1.0   1.0   1.0   1.0   0.0 |
         |  0.0   1.0   1.0   1.0   1.0  -1.0   1.0   1.0   1.0 |
         |  0.0   0.0   0.0   0.0   1.0   1.0  -1.0   1.0   1.0 |
         |  0.0   1.0   1.0   1.0   1.0   1.0   1.0  -1.0   1.0 |
         |  0.0   0.0   0.0   0.0   0.0   1.0   1.0   1.0  -1.0 |
         *                                                      *

Call Statement and Input

           N  A  NA  IDIAG  IPARM  RPARM  AUX  NAUX  BX  LDBX MBX
           |  |  |     |      |      |     |    |    |    |    |
CALL DSKFS(9, A, 33, IDIAG, IPARM, RPARM, AUX,  39 , BX , 12 , 3 )

A        =  (-1.0, -2.0, -1.0, -3.0, -2.0, -1.0, -4.0, -3.0, -2.0,
             -1.0, -4.0, -3.0, -2.0, -1.0, -5.0, -4.0, -3.0, -2.0,
             -1.0, -3.0, -2.0, -1.0, -7.0, -3.0, -5.0, -4.0, -3.0,
             -2.0, -1.0, -4.0, -3.0, -2.0, -1.0)
IDIAG    =  (1, 2, 4, 7, 11, 15, 20, 23, 30, 34)
IPARM    =  (1, 10, 0, 0, 0, . , . , . , . , 1, 0, -1, -1, -1, -1, . ,
             . , . , . , . , . , . , . , . , . )

Output

A = (-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0,
1.0, 1.0,
-1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0,
1.0, -1.0 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0,
1.0)
IDIAG =(same as input)
IPARM = (1, 10, 0, 0, 0, . , . , . , . , 1, 0, -1, -1, -1, -1, 8,
. , . , . , . , 9, 0, 0, 0, 0)
RPARM = ( . , . , . , . , . , . , . , . , . ,10^-15, . , . ,
. , . , . , 7.0, -1.0, 0.0, . , . , . , . , . , . , . )

        *                     *
        | -1.00  -2.00  -3.00 |
        | -1.00  -2.00  -3.00 |
        | -1.00  -2.00  -3.00 |
        | -1.00  -2.00  -3.00 |
        | -1.00  -2.00  -3.00 |
BX   =  | -1.00  -2.00  -3.00 |
        | -1.00  -2.00  -3.00 |
        | -1.00  -2.00  -3.00 |
        | -1.00  -2.00  -3.00 |
        |   .      .      .   |
        |   .      .      .   |
        |   .      .      .   |
        *                     *

Example 4

This example shows how to factor the first six rows and columns, referred to as matrix A1, of the 9 by 9 symmetric sparse matrix A from Example 1 and compute the determinant of A1. It uses Gaussian elimination. Input matrix A1, shown here, is stored in diagonal-out skyline storage mode. Input matrix A1 is:

                     *                              *
                     | 1.0  1.0  1.0  1.0  0.0  0.0 |
                     | 1.0  2.0  2.0  2.0  1.0  1.0 |
                     | 1.0  2.0  3.0  3.0  2.0  2.0 |
                     | 1.0  2.0  3.0  4.0  3.0  3.0 |
                     | 0.0  1.0  2.0  3.0  4.0  4.0 |
                     | 0.0  1.0  2.0  3.0  4.0  5.0 |
                     *                              *

Output matrix A1, shown here, is in LDL^T factored form with D^-1 on the diagonal, and is stored in diagonal-out skyline storage mode. Output matrix A1 is:

                     *                              *
                     | 1.0  1.0  1.0  1.0  0.0  0.0 |
                     | 1.0  1.0  1.0  1.0  1.0  1.0 |
                     | 1.0  1.0  1.0  1.0  1.0  1.0 |
                     | 1.0  1.0  1.0  1.0  1.0  1.0 |
                     | 0.0  1.0  1.0  1.0  1.0  1.0 |
                     | 0.0  1.0  1.0  1.0  1.0  1.0 |
                     *                              *

Call Statement and Input

            N   A   NA   IDIAG   IPARM   RPARM   AUX  NAUX  BX   LDBX   MBX
            |   |    |     |       |       |      |    |    |     |      |
CALL DSKFS (6 , A , 33 , IDIAG , IPARM , RPARM , AUX , 27 , BX , LDBX , MBX )

Output

Example 5

This example shows how to do a partial factorization of the 9 by 9 symmetric sparse matrix A from Example 1, where the first six rows and columns were factored in Example 4. It factors the remaining three rows and columns and computes the determinant of that part of the matrix. It uses Gaussian elimination. The input matrix, referred to as A2, shown here, is made up of the output factored matrix A1 plus the three remaining unfactored rows and columns of matrix A Matrix A2 is:

             *                                             *
             | 1.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0  0.0 |
             | 1.0  1.0  1.0  1.0  1.0  1.0  0.0  1.0  0.0 |
             | 1.0  1.0  1.0  1.0  1.0  1.0  0.0  2.0  0.0 |
             | 1.0  1.0  1.0  1.0  1.0  1.0  0.0  3.0  0.0 |
             | 0.0  1.0  1.0  1.0  1.0  1.0  1.0  4.0  0.0 |
             | 0.0  1.0  1.0  1.0  1.0  1.0  2.0  5.0  1.0 |
             | 0.0  0.0  0.0  0.0  1.0  2.0  3.0  3.0  2.0 |
             | 0.0  1.0  2.0  3.0  4.0  5.0  3.0  7.0  3.0 |
             | 0.0  0.0  0.0  0.0  0.0  1.0  2.0  3.0  4.0 |
             *                                             *

Both parts of input matrix A2 are stored in diagonal-out skyline storage mode.

Output matrix A2 is the same as output matrix A in Example 1 and is stored in diagonal-out skyline storage mode.

Call Statement and Input

            N   A   NA   IDIAG   IPARM   RPARM   AUX   NAUX  BX   LDBX   MBX
            |   |   |      |       |       |      |     |    |     |      |
CALL DSKFS (9 , A , 33 , IDIAG , IPARM , RPARM , AUX ,  27 , BX , LDBX , MBX )

Output

IDIAG =(same as output IDIAG in Example 1)
IPARM = (1, 11, 6, 0, 0, . , . , . , . , 0, . , . , . , . , . , 8,
. , . , . , . , 0, 0, 0, 0, 3)
RPARM = ( . , . , . , . , . , . , . , . , . , . , . , . , . , . ,
. , 7.0, 1.0, 0.0, . , . , . , . , . , . , . )
BX =(same as input)
LDBX =(same as input)
MBX =(same as input)

Example 6

This example shows how to solve the system Ax = b with one right-hand side for a symmetric sparse matrix A. Input matrix A, used here, is the same as factored output matrix A from Example 1, stored in profile-in skyline storage mode. It specifies Gaussian elimination, as used in Example 1. Here, output matrix A is unchanged on output and is stored in profile-in skyline storage mode.

Call Statement and Input

            N  A  NA  IDIAG  IPARM  RPARM  AUX NAUX  BX  LDBX MBX
            |  |  |     |      |      |     |   |    |    |    |
CALL DSKFS (9, A, 33, IDIAG, IPARM, RPARM, AUX, 31 , BX , 9 ,  1 )

Output

Example 7

This example shows how to factor a 9 by 9 symmetric sparse matrix A and solve the system Ax = b with four right-hand sides. It uses Cholesky decomposition. Input matrix A, shown here, is stored in profile-in skyline storage mode Matrix A is:

          *                                                    *
          | 1.0  1.0   1.0   0.0   1.0   0.0   0.0   0.0   1.0 |
          | 1.0  5.0   3.0   0.0   3.0   0.0   0.0   0.0   3.0 |
          | 1.0  3.0  11.0   3.0   5.0   3.0   3.0   0.0   5.0 |
          | 0.0  0.0   3.0  17.0   5.0   5.0   5.0   0.0   5.0 |
          | 1.0  3.0   5.0   5.0  29.0   7.0   7.0   0.0   9.0 |
          | 0.0  0.0   3.0   5.0   7.0  39.0   9.0   6.0   9.0 |
          | 0.0  0.0   3.0   5.0   7.0   9.0  53.0   8.0  11.0 |
          | 0.0  0.0   0.0   0.0   0.0   6.0   8.0  66.0  10.0 |
          | 1.0  3.0   5.0   5.0   9.0   9.0  11.0  10.0  89.0 |
          *                                                    *

Output matrix A, shown here, is in R^TR factored form with the inverse of the diagonal of R on the diagonal, and is stored in profile-in skyline storage mode. Matrix A is:

           *                                                  *
           | 1.0  1.0   1.0  0.0  1.0   0.0   0.0   0.0   1.0 |
           | 1.0   .5   1.0  0.0  1.0   0.0   0.0   0.0   1.0 |
           | 1.0  1.0  .333  1.0  1.0   1.0   1.0   0.0   1.0 |
           | 0.0  0.0   1.0  .25  1.0   1.0   1.0   0.0   1.0 |
           | 1.0  1.0   1.0  1.0   .2   1.0   1.0   0.0   1.0 |
           | 0.0  0.0   1.0  1.0  1.0  .167   1.0   1.0   1.0 |
           | 0.0  0.0   1.0  1.0  1.0   1.0  .143   1.0   1.0 |
           | 0.0  0.0   0.0  0.0  0.0   1.0   1.0  .125   1.0 |
           | 1.0  1.0   1.0  1.0  1.0   1.0   1.0   1.0  .111 |
           *                                                  *

Call Statement and Input

            N  A  NA  IDIAG  IPARM  RPARM  AUX NAUX  BX  LDBX MBX
            |  |  |     |      |      |     |   |    |    |    |
CALL DSKFS( 9, A, 34, IDIAG, IPARM, RPARM, AUX, 43 , BX , 10 , 4 )

A        =  (1.0, 1.0, 5.0, 1.0, 3.0, 11.0, 3.0, 17.0, 1.0, 3.0, 5.0,
             5.0, 29.0, 3.0, 5.0, 7.0, 39.0, 3.0, 5.0, 7.0, 9.0, 53.0,
             6.0, 8.0, 66.0, 1.0, 3.0, 5.0, 5.0, 9.0, 9.0, 11.0, 10.0,
             89.0)
IDIAG    =  (1, 3, 6, 8, 13, 17, 22, 25, 34, 35)
IPARM    =  (1, 110, 0, 1, 1, . , . , . , . , 0, . , . , . , . , . ,
             . , . , . , . , . , . , . , . , . , . )

        *                                *
        |   5.00   10.00   15.00   20.00 |
        |  15.00   30.00   45.00   60.00 |
        |  34.00   68.00  102.00  136.00 |
        |  40.00   80.00  120.00  160.00 |
BX   =  |  66.00  132.00  198.00  264.00 |
        |  78.00  156.00  234.00  312.00 |
        |  96.00  192.00  288.00  384.00 |
        |  90.00  180.00  270.00  360.00 |
        | 142.00  284.00  426.00  568.00 |
        |    .       .       .       .   |
        *                                *

Output

        *                        *
        | 1.00  2.00  3.00  4.00 |
        | 1.00  2.00  3.00  4.00 |
        | 1.00  2.00  3.00  4.00 |
        | 1.00  2.00  3.00  4.00 |
BX   =  | 1.00  2.00  3.00  4.00 |
        | 1.00  2.00  3.00  4.00 |
        | 1.00  2.00  3.00  4.00 |
        | 1.00  2.00  3.00  4.00 |
        | 1.00  2.00  3.00  4.00 |
        |  .     .     .     .   |
        *                        *