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

DSRIS--Iterative Linear System Solver for a General or Symmetric Sparse Matrix Stored by Rows

This subroutine solves a general or symmetric sparse linear system of equations, using an iterative algorithm, with or without preconditioning. The methods include conjugate gradient (CG), conjugate gradient squared (CGS), generalized minimum residual (GMRES), more smoothly converging variant of the CGS method (Bi-CGSTAB), or transpose-free quasi-minimal residual method (TFQMR). The preconditioners include an incomplete LU factorization, an incomplete Cholesky factorization (for positive definite symmetric matrices), diagonal scaling, or symmetric successive over-relaxation (SSOR) with two possible choices for the diagonal matrix: one uses the absolute values sum of the input matrix, and the other uses the diagonal obtained from the LU factorization. The sparse matrix is stored using storage-by-rows for general matrices and upper- or lower-storage-by-rows for symmetric matrices. Matrix A and vectors x and b are used:

Ax = b

where A, x, and b contain long-precision real numbers.

Syntax

On Entry

stor

indicates the form of sparse matrix A and the storage mode used, where:

If stor = 'G', A is a general sparse matrix, stored using storage-by-rows.

If stor = 'U', A is a symmetric sparse matrix, stored using upper-storage-by-rows.

If stor = 'L', A is a symmetric sparse matrix, stored using lower-storage-by-rows.

Specified as: a single character. It must be 'G', 'U', or 'L'.

init

indicates the type of computation to be performed, where:

If init = 'I', the preconditioning matrix is computed, the internal representation of the sparse matrix is generated, and the iteration procedure is performed. The coefficient matrix and preconditioner in internal format are saved in aux1.

If init = 'S', the iteration procedure is performed using the coefficient matrix and the preconditioner in internal format, stored in aux1, created in a preceding call to this subroutine with init = 'I'. You use this option to solve the same matrix for different right-hand sides, b, optimizing your performance. As long as you do not change the coefficient matrix and preconditioner in aux1, any number of calls can be made with init = 'S'.

Specified as: a single character. It must be 'I' or 'S'.

n

is the order of the linear system Ax = b and the number of rows and columns in sparse matrix A. Specified as: a fullword integer; n >= 0.

ar

is the sparse matrix A of order n, stored by rows in an array, referred to as AR. The stor argument indicates the storage variation used for storing matrix A. Specified as: a one-dimensional array, containing long-precision real numbers. The number of elements in this array can be determined by subtracting 1 from the value in IA(n+1).

ja

is the array, referred to as JA, containing the column numbers of each nonzero element in sparse matrix A. Specified as: a one-dimensional array, containing fullword integers; 1 <= (JA elements) <= n. The number of elements in this array can be determined by subtracting 1 from the value in IA(n+1).

ia

is the row pointer array, referred to as IA, containing the starting positions of each row of matrix A in array AR and one position past the end of array AR. Specified as: a one-dimensional array of (at least) length n+1, containing fullword integers; IA(i+1) >= IA(i) for i = 1, n+1.

b

is the vector b of length n, containing the right-hand side of the matrix problem. Specified as: a one-dimensional array of (at least) length n, containing long-precision real numbers.

x

is the vector x of length n, containing your initial guess of the solution of the linear system. Specified as: a one-dimensional array of (at least) length n, containing long-precision real numbers. The elements can have any value, and if no guess is available, the value can be zero.

iparm

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

• IPARM(1) controls the number of iterations.

  If IPARM(1) > 0, IPARM(1) is the maximum number of iterations allowed.

  If IPARM(1) = 0, the following default values are used:

  IPARM(1) = 300

  IPARM(2) = 4

  IPARM(4) = 4

  IPARM(5) = 1

  RPARM(1) = 10^-6

  RPARM(2) = 1

• IPARM(2) is the flag used to select the iterative procedure used in this subroutine.

  If IPARM(2) = 1, the conjugate gradient (CG) method is used. Note that this algorithm should only be used with positive definite symmetric matrices.

  If IPARM(2) = 2, the conjugate gradient squared (CGS) method is used.

  If IPARM(2) = 3, the generalized minimum residual (GMRES) method, restarted after k steps, is used.

  If IPARM(2) = 4, the more smoothly converging variant of the CGS method (Bi-CGSTAB) is used.

  If IPARM(2) = 5, the transpose-free quasi-minimal residual method (TFQMR) is used.

• IPARM(3) has the following meaning, where:

  If IPARM(2) <> 3, then IPARM(3) is not used.

  If IPARM(2) = 3, then IPARM(3) = k, where k is the number of steps after which the generalized minimum residual method is restarted. A value for k in the range of 5 to 10 is suitable for most problems.

• IPARM(4) is the flag that determines the type of preconditioning.

  If IPARM(4) = 1, the system is not preconditioned.

  If IPARM(4) = 2, the system is preconditioned by a diagonal matrix.

  If IPARM(4) = 3, the system is preconditioned by SSOR splitting with the diagonal given by the absolute values sum of the input matrix.

  If IPARM(4) = 4, the system is preconditioned by an incomplete LU factorization.

  If IPARM(4) = 5, the system is preconditioned by SSOR splitting with the diagonal given by the incomplete LU factorization.

  Note:

  The multithreaded version of DSRIS only runs on multiple threads when IPARM(4) = 1 or 2.

• IPARM(5) is the flag used to select the stopping criterion used in the computation, where the following items are used in the definitions of the stopping criteria below:
  • epsilon is the desired relative accuracy and is stored in RPARM(1).
  • x_j is the solution found at the j-th iteration.
  • r_j and r₀ are the preconditioned residuals obtained at iterations j and 0, respectively. (The residual at iteration j is given by b-Ax_j.)

  If IPARM(5) = 1, the iterative method is stopped when:

  ||r_j||₂ / ||x_j||₂ < epsilon

  Note:

  IPARM(5) = 1 is the default value assumed by ESSL if you do not specify one of the values described here; therefore, if you do not update your program to set an IPARM(5) value, you, by default, use the above stopping criterion.

  If IPARM(5) = 2, the iterative method is stopped when:

  ||r_j||₂ / ||r₀||₂ < epsilon

  If IPARM(5) = 3, the iterative method is stopped when:

  ||x_j -x_j-1||₂ / ||x_j||₂ < epsilon

  Note:

  Stopping criterion 3 performs poorly with the TFQMR method; therefore, if you specify TFQMR (IPARM(2) = 5), you should not specify stopping criterion 3.

• IPARM(6), see On Return.

Specified as: an array of (at least) length 6, containing fullword integers, where:

IPARM(1) >= 0

IPARM(2) = 1, 2, 3, 4, or 5

If IPARM(2) = 3, then IPARM(3) > 0

IPARM(4) = 1, 2, 3, 4, or 5

IPARM(5) = 1, 2, or 3 (Other values default to stopping criterion 1.)

rparm

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

RPARM(1) is the relative accuracy epsilon used in the stopping criterion. See Notes.

RPARM(2), see On Return.

RPARM(3) has the following meaning, where:

• If IPARM(4) <> 3, then RPARM(3) is not used.
• If IPARM(4) = 3, then RPARM(3) is the acceleration parameter used in SSOR. (A value in the range 0.5 to 2.0 is suitable for most problems.)

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

RPARM(1) >= 0

If IPARM(4) = 3, RPARM(3) > 0

aux1

is working storage for this subroutine, where:

If init = 'I', the working storage is computed. It can contain any values.

If init = 'S', the working storage is used in solving the linear system. It contains the coefficient matrix and preconditioner in internal format, computed in an earlier call to this subroutine.

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

naux1

is the number of doublewords in the working storage specified in aux1.

Specified as: a fullword integer, where:

In these formulas nw has the following value:

If stor = 'G', then nw = IA(n+1)-1+n.

If stor = 'U' or 'L', then nw = 2(IA(n+1)-1).

If IPARM(4) = 1, use naux1 = (3/2)nw+(7/2)n+40.

If IPARM(4) = 2, use naux1 = (3/2)nw+(9/2)n+40.

If IPARM(4) = 3, 4, or 5, then:

If IPARM(2) <> 1, use naux1 = 3nw+10n+60.

If IPARM(2) = 1, use naux1 = 3nw+(21/2)n+60.

Note:

If you receive an attention message, you have not specified sufficient auxiliary storage to achieve optimal performance, but it is enough to perform the computation. To obtain optimal performance, you need to use the amount given by the attention message.

aux2

has the following meaning:

If naux2 = 0 and error 2015 is unrecoverable, aux2 is ignored.

Otherwise, it is working storage used by this subroutine that is available for use by the calling program between calls to this subroutine.

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

naux2

is the number of doublewords in the working storage specified in aux2. Specified as: a fullword integer, where:

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

Otherwise,

If IPARM(2) = 1, use naux2 >= 4n.

If IPARM(2) = 2, use naux2 >= 7n.

If IPARM(2) = 3, use naux2 >= (k+2)n+k(k+4)+1, where k = IPARM(3).

If IPARM(2) = 4, use naux2 >= 7n.

If IPARM(2) = 5, use naux2 >= 9n.

On Return

ar

is the sparse matrix A of order n, stored by rows in an array, referred to as AR. The stor argument indicates the storage variation used for storing matrix A. The order of the elements in each row of A in AR may be changed on output.

Returned as: a one-dimensional array, containing long-precision real numbers. The number of elements in this array can be determined by subtracting 1 from the value in IA(n+1).

ja

is the array, referred to as JA, containing the column numbers of each nonzero element in sparse matrix A. These elements correspond to the arrangement of the contents of AR on output.

Returned as: a one-dimensional array, containing fullword integers; 1 <= (JA elements) <= n. The number of elements in this array can be determined by subtracting 1 from the value in IA(n+1).

x

is the vector x of length n, containing the solution of the system Ax = b. Returned as: a one-dimensional array of (at least) length n, containing long-precision real numbers.

iparm

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

IPARM(1) through IPARM(5) are unchanged.

IPARM(6) contains the number of iterations performed by this subroutine.

Returned as: a one-dimensional array of length 6, containing fullword integers.

rparm

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

RPARM(1) is unchanged.

RPARM(2) contains the estimate of the error of the solution. If the process converged, RPARM(2) <= epsilon.

RPARM(3) is unchanged.

Returned as: a one-dimensional array of length 3, containing long-precision real numbers.

aux1

is working storage for this subroutine, containing the coefficient matrix and preconditioner in internal format, ready to be passed in a subsequent invocation of this subroutine. Returned as: an area of storage, containing naux1 long-precision real numbers.

Notes

1. If you want to solve the same sparse linear system of equations multiple times using a different algorithm with the same preconditioner and using a different right-hand side each time, you get the best performance by using the following technique. Call DSRIS the first time with init = 'I'. This solves the system, and then stores the coefficient matrix and preconditioner in internal format in aux1. On the subsequent invocations of DSRIS with different right-hand sides, specify init = 'S'. This indicates to DSRIS to use the contents of aux1, saving the time to convert your coefficient matrix and preconditioner to internal format. If you use this technique, you should not modify the contents of aux1 between calls to DSRIS.

   In some cases, you can specify a different algorithm in IPARM(2) when making calls with init = 'S'. (See Example 2.) However, DSRIS sometimes needs different information in aux1 for different algorithms. When this occurs, DSRIS issues an attention message, continues processing the computation, and then resets the contents of aux1. Your performance is not improved in this case, which is functionally equivalent to calling DSRIS with init = 'I'.

2. If you use the CG method with init = 'I', you must use the CG method when you specify init = 'S'. However, if you use a different method with init = 'I', you can use any other method, except CG, when you specify init = 'S'.
3. These subroutines accept lowercase letters for the stor and init arguments.
4. Matrix A, vector x, and vector b must have no common elements; otherwise, results are unpredictable.
5. In this subroutine, a value of RPARM(1) = 0 is permitted to force the solver to evaluate exactly IPARM(1) iterations. The algorithm computes a sequence of approximate solution vectors x that converge to the solution. The iterative procedure is stopped when the selected stopping criterion is satisfied or when more than the maximum number of iterations (in IPARM(1)) is reached.

   For the stopping criteria specified in IPARM(5), the relative accuracy epsilon (in RPARM(1)) must be specified reasonably (10^-4 to 10^-8). If you specify a larger epsilon, the algorithm takes fewer iterations to converge to a solution. If you specify a smaller epsilon, the algorithm requires more iterations and computer time, but converges to a more precise solution. If the value you specify is unreasonably small, the algorithm may fail to converge within the number of iterations it is allowed to perform.

6. For a description of how sparse matrices are stored by rows, see Storage-by-Rows.
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

The linear system:

Ax = b

is solved using one of the following methods: conjugate gradient (CG), conjugate gradient squared (CGS), generalized minimum residual (GMRES), more smoothly converging variant of the CGS method (Bi-CGSTAB), or transpose-free quasi-minimal residual method (TFQMR), where:

A is a sparse matrix of order n. The matrix is stored in arrays AR, IA, and JA. If it is general, it is stored by rows. If it is symmetric, it can be stored using upper- or lower-storage-by-rows.

x is a vector of length n.

b is a vector of length n.

One of the following preconditioners is used:

• an incomplete LU factorization
• an incomplete Cholesky factorization (for positive definite symmetric matrices)
• diagonal scaling
• symmetric successive over-relaxation (SSOR) with two possible choices for the diagonal matrix:
  • the absolute values sum of the input matrix
  • the diagonal obtained from the LU factorization

See references [36], [56], [82], [86], [89], and [95].

When you call this subroutine to solve a system for the first time, you specify init = 'I'. After that, you can solve the same system any number of times by calling this subroutine each time with init = 'S'. These subsequent calls use the coefficient matrix and preconditioner, stored in internal format in aux1. You optimize performance by doing this, because certain portions of the computation have already been performed.

Error Conditions

Resource Errors

Error 2015 is unrecoverable, naux2 = 0, and unable to allocate work area.

Computational Errors

The following errors, with their corresponding return codes, can occur in this subroutine. For details on error handling, see What Can You Do about ESSL Computational Errors?.

• For error 2110, return code 1 indicates that the subroutine exceeded IPARM(1) iterations without converging. Vector x contains the approximate solution computed at the last iteration.
• For error 2130, return code 2 indicates that the incomplete LU factorization of A could not be completed, because one pivot was 0.
• For error 2124, the subroutine has been called with init = 'S', but the data contained in aux1 was computed for a different algorithm. An attention message is issued. Processing continues, and the contents of aux1 are reset correctly.
• For error 2134, return code 3 indicates that the data contained in aux1 is not consistent with the input sparse matrix. The subroutine has been called with init = 'S', and aux1 contains an incomplete factorization and internal data storage for the input matrix A that was computed by a previous call to the subroutine when init = 'I'. This error indicates that aux1 has been modified since the last call to the subroutine, or that the input matrix is not the same as the one that was factored. If the default action has been overridden, the subroutine can be called again with the same parameters, with the exception of IPARM(4) = 1 or 4.
• For error 2131, return code 4 indicates that the matrix is singular, because all elements in one row of the matrix contain zero.
• For error 2129, return code 5 indicates that the matrix is not positive definite.
• For error 2128, return code 8 indicates an internal ESSL error. Please contact your IBM Representative.

Input-Argument Errors

1. n < 0
2. stor <> 'G', 'U', or 'L'
3. init <> 'I' or 'S'
4. IA(n+1) < 1
5. IA(i+1)-IA(i) < 0, for any i = 1, n
6. IPARM(1) < 0
7. IPARM(2) <> 1, 2, 3, 4, or 5
8. IPARM(3) <= 0 and IPARM(2) = 3
9. IPARM(4) <> 1, 2, 3, 4, or 5
10. RPARM(1) < 0
11. RPARM(3) <= 0 and IPARM(4) = 3
12. naux1 is too small--that is, less than the minimum required value. Return code 6 is returned if error 2015 is recoverable.
13. Error 2015 is recoverable or naux2<>0, and naux2 is too small--that is, less than the minimum required value. Return code 7 is returned for naux2 if error 2015 is recoverable.

Example 1

This example finds the solution of the linear system Ax = b for the sparse matrix A, which is stored by rows in arrays AR, IA, and JA. The system is solved using the Bi-CGSTAB algorithm. The iteration is stopped when the norm of the residual is less than the given threshold specified in RPARM(1). The algorithm is allowed to perform 20 iterations. The process converges after 9 iterations. Matrix A is:

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

Call Statement and Input

            STOR INIT   N   AR   JA   IA   B   X   IPARM   RPARM   AUX1   NAUX1   AUX2   NAUX2
             |     |    |   |    |    |    |   |     |       |      |       |      |       |
CALL DSRIS( 'G' , 'I' , 9 , AR , JA , IA , B , X , IPARM , RPARM , AUX1 ,   98  , AUX2 ,   63 )

AR       =  (2.0, 2.0, -1.0, 1.0, 2.0, 1.0, 2.0, -1.0, 1.0, 2.0, -1.0,
             1.0, 2.0, -1.0, 1.0, 2.0, -1.0, 1.0, 2.0, -1.0, 1.0, 2.0)
JA       =  (1, 2, 3, 2, 3, 1, 4, 5, 4, 5, 6, 5, 6, 7, 6, 7, 8, 7, 8,
             9, 8, 9)
IA       =  (1, 2, 4, 6, 9, 12, 15, 18, 21, 23)
B        =  (2.0, 1.0, 3.0, 2.0, 2.0, 2.0, 2.0, 2.0, 3.0)
X        =  (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
IPARM(1) =  20
IPARM(2) =  4
IPARM(3) =  0
IPARM(4) =  1
IPARM(5) =  10
RPARM(1) =  1.D-7
RPARM(3) =  1.0

Output

X        =  (1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0)
IPARM(6) =  9
RPARM(2) =  0.29D-16

Example 2

This example finds the solution of the linear system Ax = b for the same sparse matrix A used in Example 1. It also uses the same right-hand side in b and the same initial guesses in x. However, the system is solved using a different algorithm, conjugate gradient squared (CGS). Because INIT is 'S', the best performance is achieved. The iteration is stopped when the norm of the residual is less than the given threshold specified in RPARM(1). The algorithm is allowed to perform 20 iterations. The process converges after 9 iterations.

Call Statement and Input

            STOR INIT   N   AR   JA   IA   B   X   IPARM   RPARM   AUX1   NAUX1   AUX2   NAUX2
             |     |    |   |    |    |    |   |     |       |      |       |      |       |
CALL DSRIS( 'G' , 'S' , 9 , AR , JA , IA , B , X , IPARM , RPARM , AUX1 ,   98  , AUX2 ,   63 )

Output

X        =  (1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0)
IPARM(6) =  9
RPARM(2) =  0.42D-19

Example 3

This example finds the solution of the linear system Ax = b for the sparse matrix A, which is stored by rows in arrays AR, IA, and JA. The system is solved using the two-term conjugate gradient method (CG), preconditioned by incomplete LU factorization. The iteration is stopped when the norm of the residual is less than the given threshold specified in RPARM(1). The algorithm is allowed to perform 20 iterations. The process converges after 1 iteration. Matrix A is:

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

Call Statement Input

            STOR INIT   N   AR   JA   IA   B   X   IPARM   RPARM   AUX1   NAUX1   AUX2   NAUX2
             |     |    |   |    |    |    |   |     |       |      |       |      |       |
CALL DSRIS( 'G' , 'I' , 9 , AR , JA , IA , B , X , IPARM , RPARM , AUX1 ,  223  , AUX2 ,   36 )

AR       =  (2.0, -1.0, 2.0, -1.0, -1.0, 2.0, -1.0, -1.0, 2.0, -1.0,
             -1.0, 2.0, -1.0, -1.0, 2.0, -1.0, -1.0, 2.0, -1.0, -1.0,
             2.0, -1.0, 2.0)
JA       =  (1, 3, 2, 4, 1, 3, 5, 2, 4, 6, 3, 5, 7, 4, 6, 8, 5, 7, 9,
             6, 8, 7, 9)
IA       =  (1, 3, 5, 8, 11, 14, 17, 20, 22, 24)
B        =  (1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0)
X        =  (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
IPARM(1) =  20
IPARM(2) =  1
IPARM(3) =  0
IPARM(4) =  4
IPARM(5) =  1
RPARM(1) =  1.D-7
RPARM(3) =  1.0

Output

X        =  (1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0)
IPARM(6) =  1
RPARM(2) =  0.16D-15

Example 4

This example finds the solution of the linear system Ax = b for the same sparse matrix A used in Example 3. However, matrix A is stored using upper-storage-by-rows in arrays AR, IA, and JA. The system is solved using the generalized minimum residual (GMRES), restarted after 5 steps and preconditioned with SSOR splitting. The iteration is stopped when the norm of the residual is less than the given threshold specified in RPARM(1). The algorithm is allowed to perform 20 iterations. The process converges after 12 iterations.

Call Statement Input

            STOR INIT   N   AR   JA   IA   B   X   IPARM   RPARM   AUX1   NAUX1   AUX2   NAUX2
             |     |    |   |    |    |    |   |     |       |      |       |      |       |
CALL DSRIS( 'U' , 'I' , 9 , AR , JA , IA , B , X , IPARM , RPARM , AUX1 ,  219  , AUX2 ,  109 )

AR       =  (2.0, -1.0, 2.0, -1.0, 2.0, -1.0, 2.0, -1.0, 2.0, -1.0,
             2.0, -1.0, 2.0, -1.0, 2.0, 2.0)
JA       =  (1, 3, 2, 4, 3, 5, 4, 6, 5, 7, 6, 8, 7, 9, 8, 9)
IA       =  (1, 3, 5, 7, 9, 11, 13, 15, 16, 17)
B        =  (1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0)
X        =  (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
IPARM(1) =  20
IPARM(2) =  3
IPARM(3) =  5
IPARM(4) =  3
IPARM(5) =  1
RPARM(1) =  1.D-7
RPARM(3) =  2.0

Output

X        =  (1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0)
IPARM(6) =  12
RPARM(2) =  0.33D-7