The Y12M documentation: C ================================= C Documentation of subroutine Y12MA C ================================= C C 1. Purpose C C C Y12MA solves sparse systems of linear algebraic C equations by Gaussian elimination. The subroutine is a C "black box subroutine" designed to solve efficiently C problems which contain only one system with a single C right hand side. The number of the input parameters is C minimized. The user must assign values only to NN, NN1, C N, Z, A, SNR, RNR, IHA and B according to the rules C described in Section 2.4 (see below). It is extremely C easy to modify the subroutine to the cases: (a) a C sequence of systems with the same matrix is to be C solved (note that one system with many right hand sides C can be rewritten as a sequence of systems with the same C matrix), (b) a sequence of systems whose matrices are C different but of the same structure is to be solved and C (c) a sequence of systems whose matrices are of the C same structure and some of them are the same is to be C solved. These cases are defined as case (ii), case C (iii) and case (iv) in Section 1.1. "Scope of the C Y12M". The recommendations in Section 1.3.6 should be C followed in order to modify the subroutine Y12MA for C the above cases. If a sequence of systems whose C matrices are of different structure (this case appears C as case (v) in Section 1.1) is to be solved then C subroutine Y12MA should be called in the solution of C each system in the sequence. C C 2. Calling sequence and declaration of the parameters C C The subroutine is written in FORTRAN and it has been C extensively tested with the FOR and FTN compilers on a C UNIVAC 1100/82 computer, at the Regional Computing C Centre at the University of Copenhagen (RECKU). Many C examples have been run on an IBM 3033 computer at the C Northern Europe University Computing Centre (NEUCC) and C on a CDC Cyber 173 computer at the Regional Computing C Centre at the University of Aarhus (RECAU). Two C different versions are available: a single precision C version named Y12MAE and a double precision version C named Y12MAF. The calls of these two versions and the C declarations of the parameters are as follows. C C A). Single precision version: Y12MAE C C SUBROUTINE Y12MAE(N, Z, A, SNR, NN, RNR, NN1, C 1 PIVOT, HA, IHA, AFLAG, IFLAG, B, IFAIL) C C REAL A(NN), PIVOT(N), B(N), AFLAG(8) C INTEGER SNR(NN), RNR(NN1), HA(IHA,11), IFLAG(10) C INTEGER N, Z, NN, NN1, IHA, IFAIL C C B). Double precision version: Y12MAF C C SUBROUTINE Y12MAF(N, Z, A, SNR, NN, RNR, NN1, C 1 PIVOT, HA, IHA, AFLAG, IFLAG, B, IFAIL) C C DOUBLE PRECISION A(NN), PIVOT(N), B(N), AFLAG(8) C INTEGER SNR(NN), RNR(NN1), HA(IHA,11), IFLAG(10) C INTEGER N, Z, NN, NN1, IHA, IFAIL C C These two versions can be used on many other computers C also. However some alterations may be needed and/or may C ensure greater efficiency of the performance of the C subroutine. For example, it will be much more efficient C to declare arrays SNR, RNR and (if possible) HA as C INTEGER*2 arrays on some IBM installations. C C C 3. Method C C The system Ax = b is solved by Gaussian elimination. C Pivotal interchanges are used in an attempt to preserve C both the stability of the computations and the sparsity C of the original matrix. In this way a decomposition LU C = PAQ is normally calculated. P and Q are permutation C matrices, L is a lower triangular matrix, U is an upper C triangular matrix. The right hand side vector b is also C modified during the decomposition so that the vector C c=L**(-1)*P*b is available after the decomposition C stage. In this way there is no need to store the C non-zero elements of matrix L. Therefore these elements C are not stored. This leads to a reduction of the C storage requirements (the length of arrays A and SNR C can be decreased). An approximation to the solution C vector is found by solving U*tr(T)*x=c, where tr(T) is C the transpose. The subroutine calculates the rate at C which the elements of the matrix grow during the C decomposition (see parameter AFLAG(5) below) and the C minimal in absolute value pivotal element (see C parameter AFLAG(8) below). These two numbers can be C used to decide whether the computed approximation is C acceptable or not. Positive values of a special C parameter IFAIL indicate that the subroutine is unable C to solve the problem. The error diagnostics, given in C Section 7, describe the probable cause. C C 4. Parameters of the subroutine C C C N - INTEGER On entry N must contain the number of C equations in the system Ax=b. Unchanged on C exit. C C Z - INTEGER. On entry Z must contain the number of C non-zero elements in the coefficient matrix A C of the system Ax = b. Unchanged on exit. C C A - REAL (in the single precision version Y12MAE) C or DOUBLE PRECISION (in the double precision C version Y12MAF) array of length NN (see below). C On entry the first Z locations of array A must C contain the non-zero elements of the C coefficient matrix A of the system Ax = b. The C order of the non-zero elements may be C completely arbitrary. The content of array A is C modified by subroutine Y12MA. On successful C exit array A will contain the non-zero elements C of the upper triangular matrix U (without the C diagonal elements of matrix U which can be C found in array PIVOT, see below). C C SNR - INTEGER array of length NN (see below). On C entry SNR(j), j = 1(1)Z, must contain the C column number of the non-zero element stored in C A(j). The content of array SNR is modified by C subroutine Y12MA. On successful exit array SNR C will contain the column numbers of the non-zero C elements of the upper triangular matrix U C (without the column numbers of the diagonal C elements of matrix U). C C NN - INTEGER. On entry NN must contain the length of C arrays A and SNR. Restriction: NN.GE.2*Z. C Recommended value: 2*Z.LE.NN.LE.3*Z. Unchanged C on exit. C C RNR - INTEGER array of length NN1 (see below). On C entry RNR(i), i = 1(1)Z, must contain the row C number of the non-zero element stored in A(i). C The content of array RNR is modified by C subroutine Y12MA. On successful exit all C components of array RNR will normally be zero. C C NN1 - INTEGER. On entry NN1 must contain the length C of the array RNR. Restriction: NN1.ge.Z. C Recommended value: 2*Z.LE.NN1.LE.3*Z. C Unchanged on exit. C C PIVOT - REAL (in the single precision version Y12MAE) C or DOUBLE PRECISION (in the double precision C version Y12MAF) array of length N. The content C of array PIVOT is modified by subroutine Y12MA. C On successful exit array PIVOT will contain the C pivotal elements (the diagonal elements of C matrix U). This means that a small element (or C small elements) in array PIVOT on exit may C indicate numerical singularity of the C coefficient matrix A. Note that the smallest in C absolute value element in array PIVOT is also C stored in AFLAG(8), see below. C C HA - INTEGER two-dimensional array. The length of C the first dimension is IHA (see below). The C length of the second dimension is 11. This C array is used as a work space by subroutine C Y12MA. C C IHA - INTEGER. On entry IHA must contain the length C of the first dimension of array HA. C Restriction: IHA.GE.N. Unchanged on exit. C C AFLAG - REAL (in the single precision version Y12MAE) C or DOUBLE PRECISION (in the double precision C version Y12MAF) array of length 8. The content C of array AFLAG is modified by subroutine Y12MA. C The content of the components of this array can C be described as follows. C C AFLAG(1) - Stability factor. An element can be C chosen as pivotal element only if C this element is larger (in absolute C value) than the absolute value of C the largest element in its row C divided by AFLAG(1). Subroutine C Y12MA sets AFLAG(1) = 16. This C value has been found satisfactory C for a wide range of test-examples. C C AFLAG(2) - Drop-tolerance. An element, which in C the process of the computations C becomes smaller (in absolute value) C than the drop-tolerance, is removed C from array A (and its column and row C numbers are removed from arrays SNR C and RNR). Subroutine Y12MA sets C AFLAG(2) = 1.0E-12. This value has C been found satisfactory for a wide C range of test-matrices. By this C choice it is assumed that the matrix C is not too badly scaled and that the C magnitude of the elements is 1. C C AFLAG(3) - The subroutine will stop the C computations when the growth factor C (parameter AFLAG(5), see below) C becomes larger than AFLAG(3). Our C experiments show that if AFLAG(3) > C 1.0E6 then the solution is normally C quite wrong. Therefore subroutine C Y12MA sets AFLAG(3) = 1.0E6. C C AFLAG(4) - The subroutine will stop the C computations when the absolute value C of a current pivotal element is C smaller than AFLAG(4)*AFLAG(6) C (where the absolute value of the C largest (in absolute value) element C of the original matrix is stored in C AFLAG(6), see below). Our C experiments show that C AFLAG(4)=1.0E-12 will normally be a C good choice. Therefore subroutine C Y12MA sets AFLAG(4)=1.0E-12. C C AFLAG(5) - Growth factor. After each stage of C the elimination subroutine Y12MA C sets AFLAG(5) = AFLAG(7)/AFLAG(6) C (the description of parameters C AFLAG(6) and AFLAG(7) is given C below). Large values of AFLAG(5) C indicate that an appreciable error C in the computed solution is C possible. In an extreme case when C AFLAG(5) becomes larger than 1.0E6, C see parameter AFLAG(3), subroutine C Y12MA will stop the computations in C an attempt to prevent overflows. C C AFLAG(6) - Subroutine Y12MA sets AFLAG(6) equal C to max(abs(a(i,j))), i = 1(1)N, j = C 1(1)N. C C AFLAG(7) - Subroutine Y12MA sets AFLAG(7) equal C to max(abs(a(i,j;s))), 1.lt.s.lt.k, C after each step k, k=1(1)N-1, of C the elimination. C C AFLAG(8) - Subroutine Y12MA sets AFLAG(8) equal C to min(abs(a(i,i))), i= 1(1)N. This C means that the minimal pivotal C element (in absolute value) will be C stored in AFLAG(8) on successful C exit. Small values of AFLAG(8) C indicate numerical singularity of C the original matrix. We advise the C user to check this parameter on exit C very carefully. C C IFLAG - INTEGER array of length 10. The content of this C array is modified by subroutine Y12MA. The C content of the components of this array can be C described as follows. C C IFLAG(1) - Subroutine Y12MA uses IFLAG(1) as a C work space. C C IFLAG(2) - Subroutine Y12MA sets IFLAG(2) = 3. C This means that at each stage of the C Gaussian elimination (except the C last one) the pivotal search is C carried out in the 3 rows which have C least numbers of non-zero elements. C C IFLAG(3) - Subroutine Y12MA sets IFLAG(3) = 1. C This means that the pivotal strategy C for general matrices will be used. C For some special matrices (for C example positive definite) this will C be inefficient. Subroutine Y12MA C can easily be modified for such C matrices; only the statement C IFLAG(3) = 1 should be changed (for C positive definite matrices e.g. C IFLAG(3) = 2 can be used). More C details about the use of this C parameter with special matrices are C given in Section 1.3.1 C C IFLAG(4) - Subroutine Y12MA sets IFLAG(4) = 0. C This is the best choice in the case C where only one system Ax = b is to C be solved, see more details in C Section 1.3. C C IFLAG(5) - Subroutine Y12MA sets IFLAG(5) = 1. C This means that the non-zero C elements of the lower triangular C matrix L will be removed during the C decomposition stage, see more C details in Section 1.3.3. C C IFLAG(6) - On successful exit IFLAG(6) will be C equal to the number of "garbage" C collections in the row ordered list. C If IFLAG(6) is large then it is C better to choose a larger value of C NN in the next calls of subroutine C Y12MA with the same or a similar C matrix A. This will lead to a C reduction in the computing time. C C IFLAG(7) - On successful exit IFLAG(7) will be C equal to the number of "garbage" C collections in the column ordered C list. If IFLAG(7) is large then it C is better to choose a larger value C of NN1 in the next calls of C subroutine Y12MA with the same or a C similar matrix A. This will lead to C a reduction in the computing time. C C IFLAG(8) - On successful exit IFLAG(8) will be C equal to the maximal number of C non-zero elements kept in array A at C any stage of the Gaussian C elimination. If IFLAG(8) is much C smaller then NN (or NN1) then the C length NN (or NN1) can be chosen C smaller in the next calls of C subroutine Y12MA with the same or a C similar matrix A. This will lead to C a reduction of the storage needed. C C IFLAG(9) - This parameter is ignored by C subroutine Y12MA. It will be used C in some of the other subroutines C when IFLAG(4) = 1 is specified on C entry. C C IFLAG(10)- This parameter is ignored by C subroutine Y12MA. It will be used C in some of the other subroutines C when IFLAG(4) = 1 is specified on C entry. C C B - REAL (in the single precision version Y12MAE) C or DOUBLE PRECISION (in the double precision C version Y12MAF) array of length N. On entry the C right hand side vector b of the system Ax = b C must be stored in array B. The content of array C B is modified by subroutine Y12MA. On C successful exit the computed solution vector C will be stored in array B. C C IFAIL - Error diagnostic parameter. The content of C parameter IFAIL is modified by subroutine C Y12MA. On exit IFAIL = 0 if the subroutine has C not detected any error. Positive values of C IFAIL on exit show that some error has been C detected by the subroutine. Many of the error C diagnostics are common for all subroutines in C the package. Therefore the error diagnostics C are listed in a separate section, Section 7, of C this book. We advise the user to check the C value of this parameter on exit. C C 5. Error diagnostics C C Error diagnostics are given by positive values of the C parameter IFAIL (see above). We advise the user to C check carefully the value of this parameter on exit. C The error messages are listed in Section 7. C C 6. Auxiliary subroutines C C Y12MA calls three other subroutines: Y12MB, Y12MC and C Y12MD. C C 7. Timing C C The time taken depends on the order of the matrix C (parameter N), the number of the non-zero elements in C the matrix (parameter Z), the magnitude of the non-zero C elements and their distribution in the matrix. C C 8. Storage C C There are no internally declared arrays. C C 9. Accuracy C C It is difficult to evaluate the accuracy of the C computed solution. Large values of parameter AFLAG(5), C the growth factor, indicate unstable computations C during the decomposition stage. Small values of C parameter AFLAG(8) can be considered as a signal for C numerical singularity. We must emphasize here that C normally much more reliable evaluations of the accuracy C achieved can be found by the use of iterative C refinement, i.e. by the use of subroutine Y12MF. By the C use of the latter subroutine the computing time will C often be reduced too. However, the storage requirements C may be increased sometimes. C C 10. Some remarks C C Remark 1 Subroutine Y12MA may also be considered as C an example of how to initialize the first 4 C components of array AFLAG and the first 5 C components of array IFLAG when only one C system with a single right hand side is to C be solved. An efficient use of the C subroutines Y12MB, Y12MC and Y12MD for the C other cases discussed in Section 1.1 may be C achieved by a single modification of Y12MA C according to the rules given in Section C 1.3.6 (a change of only one or two C statements in Y12MA is needed to obtain such C a modification). C C Remark 2 The last 4 components of array AFLAG and the C last 5 components of array IFLAG can be used C as output parameters. The user is not C obliged to do this. However, we recommend C the check of these parameters on exit. C ================================= C Documentation of subroutine Y12MB C ================================= C C 1. Purpose C C Y12MB prepares a system of linear algebraic equations C to be factorized (by subroutine Y12MC) and solved (by C subroutine Y12MD). Each call of subroutine Y12MC must C be preceded by a call of Y12MB. It is very convenient C to perform some operations on the non-zero elements of C the coefficient matrix between calls of Y12MB and C Y12MC. For example, when iterative refinement is C carried out the information needed in the calculation C of the residual vector is copied in some extra arrays C between the calls of Y12MB and Y12MC. It is also very C easy to perform any kind of scaling after the call of C Y12MB. C C 2. Calling sequence and declaration of the parameters C C The subroutine is written in FORTRAN and has been C extensively tested with the FOR and FTN compilers on a C UNIVAC 1100/82 computer at the Regional Computing C Centre at the University of Copenhagen (RECKU). Many C examples have been run on an IBM 3033 computer at the C Northern Europe University Computing Centre (NEUCC) and C on a CDC Cyber 173 computer at the Regional Computing C Centre at the University of Aarhus (RECAU). Two C different versions are available: a single precision C version Y12MBE and a double precision version Y12MBF. C The calls of these two versions and the declaration of C the parameters are as follows: C C A) Single precision version Y12MBE C C SUBROUTINE Y12MBE(N, Z, A, SNR, NN, RNR, NN1, C 1 HA, IHA, AFLAG, IFLAG, IFAIL) C C REAL A(NN), AFLAG(8) C INTEGER SNR(NN), RNR(NN1), HA(N,11), IFLAG(10) C INTEGER N, Z, IHA, IFAIL C C B) Double precision version Y12MBF C C SUBROUTINE Y12MBF(N, Z, A, SNR, NN, RNR, NN1, C 1 HA, IHA, AFLAG, IFLAG, IFAIL) C C DOUBLE PRECISION A(NN), AFLAG(8) C INTEGER SNR(NN), RNR(NN1), HA(N,11), IFLAG(10) C INTEGER N, Z, IHA, IFAIL C C These two versions can be used on many other computers C also. However some alterations may be needed and/or may C ensure greater efficiency of the performance of the C subroutine. For example, it will be much more efficient C to declare arrays SNR, RNR and (if possible) HA as C INTEGER*2 arrays on some IBM installations. C C C 3. Method C C The non-zero elements of matrix A are ordered by rows C and stored in the first Z positions of array A. The C order of the non-zero elements within a row is C arbitrary. The column numbers of the non-zero elements C are stored in the first Z positions of array SNR so C that if A(J)=a(i,j), J=1(1)Z, then SNR(J)=j. The row C numbers of the non-zero elements are stored in the C first Z positions of array RNR so that the row numbers C of the non-zero elements in the first column of matrix C A are located before the row numbers of the non-zero C elements in the second column of matrix A, the row C numbers of the non-zero elements in the second column C of matrix A are located before the row numbers of the C non-zero elements in the third column of matrix A and C so on. Some additional information, e.g. about the row C starts, row ends, column starts and column ends, is C stored in the work array HA. This storage scheme has C been proposed by Gustavson. C C 4. Parameters of the subroutine C C C N - INTEGER. On entry N must contain the number of C equations in the system Ax=b. Unchanged on C exit. C C Z - INTEGER. On entry Z must contain the number of C non-zero elements in the coefficient matrix A C of the system Ax = b. Unchanged on exit. C C A - REAL (in the single precision version Y12MBE) C or DOUBLE PRECISION (in the double precision C version Y12MBF) array of length NN (see below). C On entry the first Z locations of array A must C contain the non-zero elements of the C coefficient matrix A of the system Ax = b. The C order of the non-zero elements may be C completely arbitrary. The content of array A is C modified by subroutine Y12MB. On successful C exit the first Z positions of array A will C contain the non-zero elements of the matrix A C ordered by rows. C C SNR - INTEGER array of length NN (see below). On C entry SNR(j), j = 1(1)Z, must contain the C column number of the non-zero element stored in C A(j). The content of array SNR is modified by C subroutine Y12MB. On successful exit SNR(j), C j=1(1)Z, will contain the column number of the C non-zero element stored in A(j). C C NN - INTEGER. On entry NN must contain the length of C arrays A and SNR. Restriction: NN.ge.2*Z. C Recommended value: 2*Z.le.NN.le.3*Z. See also C the description of NN in Y12MC. Unchanged on C exit. C C RNR - INTEGER array of length NN1 (see below). On C entry RNR(i), i = 1(1)Z, must contain the row C number of the non-zero element stored in A(i). C The content of array RNR is modified by C subroutine Y12MB. On successful exit the row C numbers of the non-zero elements of matrix A C are stored in the first Z positions of array C RNR, so that the row numbers of the non-zero C elements in the first column of matrix A are C before the row numbers of the non-zero elements C in the second column of matrix A, the row C numbers of the non-zero elements in the second C column of matrix A are before the row numbers C of the non-zero elements in the third column of C matrix A and so on. This means that on exit the C row number of the non-zero element stored in C A(i), i = 1(1)Z, is not stored in RNR(i), in C general. C C NN1 - INTEGER. On entry NN1 must contain the length C of the array RNR. Restriction: NN1.ge.Z. C Recommended value: 2*Z.le.NN1.le.3*Z. See also C the desription of NN1 in Y12MC. Unchanged on C exit. C C HA - INTEGER two-dimensional array. The length of C the first dimension of HA is IHA (see below). C The length of the second dimension of HA is 11. C The contents of some elements of this array are C modified by subroutine Y12MB, the contents of C the others are ignored. On successful exit: C C (1) HA(i,1) contains the position in array A C where the first element of row i, i C = 1(1)N, is stored. C (2) HA(i,3) contains the position in array A C where the last element of row i, i C = 1(1)N, is stored (all non-zero C elements of row i are located C between HA(i,1) and HA(i,3) C compactly; the number of non-zero C elements in row i is HA(i,3) - C HA(i,1) + 1). C (3) HA(i,4) contains the position in array RNR C where the row number of the first C non-zero element of column j is C stored (j = 1(1)N). C (4) HA(j,6) contains the position in array RNR C where the row number of the last C element in column j, j = 1(1)N, is C stored (all row numbers of the C non-zero elements of column j are C located between HA(j,4) and HA(j,6) C compactly, the number of non-zero C elements in column j is C HA(j,6)-HA(j,4)+1). C C (5) Some information needed in the C pivotal search during the C decomposition is stored in columns C 7,8 and 11 of array HA. C C (6) The other columns of HA are either C ignored by Y12MB or used as a work C space. C C IHA - INTEGER. On entry IHA must contain the length C of the first dimension of array HA. C Restriction: IHA.ge.N. Unchanged on exit. C C AFLAG - REAL (in the single precision version Y12MBE) C or DOUBLE PRECISION (in the double precision C version Y12MBF) array of length 8. The contents C of all locations of array AFLAG, except C AFLAG(6), are ignored by subroutine Y12MB. The C content of AFLAG(6) is modified by Y12MB. On C successful exit AFLAG(6) contains C max(abs(i,j)). C C IFLAG - INTEGER array of length 10. The user should set C IFLAG(1).ge.0 before the call of package Y12M C (i.e. before the first call of a subroutine of C this package). IFLAG(1) is used in the error C checks by the subroutines and should not be C altered by the user between any two successive C calls of subroutines of the package. IFLAG(1) C will be equal to -1 on successful exit from C subroutine Y12MB. Thus, the content of IFLAG(1) C is modified by subroutine Y12MB. On entry C IFLAG(4) must contain 0, 1 or 2. IFLAG(4) = 0 C is the best choice when only one system is to C be solved, when the first system of a sequence C of systems with the same matrix is to be solved C and when any system of a sequence of systems C whose matrices are of different structure is to C be solved. IFLAG(4) = 1 is the best choice when C the first system in a sequence of systems with C the same structure is to be solved. IFLAG(4) = C 2 is the best choice when any system after the C first one in a sequence of systems whose C matrices are of the same structure is to be C solved. The content of IFLAG(4) is unchanged by C subroutine Y12MB. The other locations of IFLAG C are either ignored by Y12MB or used as a work C space. C C IFAIL - Error diagnostic parameter. The content of C parameter IFAIL is modified by subroutine C Y12MB. On exit IFAIL = 0 if the subroutine has C not detected any error. Positive values of C IFAIL on exit show that some error has been C detected by the subroutine. Many of the error C diagnostics are common for all subroutines in C the package. Therefore the error diagnostics C are listed in a separate section, Section 7, of C this book. We advise the user to check the C value of this parameter on exit. C C 5. Error diagnostics C C Error diagnostics are given by positive values of the C parameter IFAIL (see above). We advise the user to C check carefully the value of this parameter on exit. C The error messages are listed in Section 7. C C 6. Auxiliary subroutines C C None. C C 7. Timing C C The time taken depends on the number of non-zero C elements (parameter Z) in matrix A. C C 8. Storage C C There are no internally declared arrays. C C 9. Some remarks C C Remark 1 The use of subroutine Y12MB is followed by the C use of Y12MC and Y12MD. Subroutines Y12MA and C Y12MF can be considered as examples of how to C use Y12MB, Y12MC and Y12MD. C C Remark 2 The contents of N, Z, A, SNR, NN, RNR, NN1, C columns 1, 3, 4, 6, 7, 8 and 11 of HA, IHA, C AFLAG(6), IFLAG(1), IFLAG(4) and IFAIL should C not be altered between calls of Y12MB and C Y12MC. C C Remark 3 If IFAIL > 0 on exit of Y12MB there is no C sense in calling Y12MC and Y12MD. Therefore C the computations should be stopped in this C case and one should investigate (using the C error diagnostics given in Section 7) why C Y12MB has assigned a positive value to IFAIL. C ================================= C Documentation of subroutine Y12MC C ================================= C C 1. Purpose C C Y12MC decomposes a matrix A into two triangular C matrices L and U. Each call of subroutine Y12MC should C be preceded by a call of Y12MB. Subroutine Y12MD is C normally called one or several times after the call of C Y12MC. C C 2. Calling sequence and declaration of the parameters C C The subroutine is written in FORTRAN and has been C extensively tested with the FOR and FTN compilers on C the UNIVAC 1100/82 computer at the Regional Computing C Centre at the University of Copenhagen (RECKU). Many C examples have been run on an IBM 3033 computer at the C Northern Europe University Computing Centre (NEUCC) and C on a CDC Cyber 173 computer at the Regional Computing C Centre at the University of Aarhus (RECAU). Two C different versions are available: a single precision C version named Y12MCE and a double precision version C named Y12MCF. The calls of these two versions and the C declaration of the parameters are as follows. C C A). Single precision version: Y12MCE C C SUBROUTINE Y12MCE(N, Z, A, SNR, NN, RNR, NN1, C 1 PIVOT, B, HA, IHA, AFLAG, IFLAG, IFAIL) C C REAL A(NN), PIVOT(N), B(N), AFLAG(8) C INTEGER SNR(NN), RNR(NN1), HA(IHA,11), IFLAG(10) C INTEGER N, Z, NN, NN1, IHA, IFAIL C C B). Double precision version: Y12MCF C C SUBROUTINE Y12MCF(N, Z, A, SNR, NN, RNR, NN1, C 1 PIVOT, B, HA, IHA, AFLAG, IFLAG, IFAIL) C C DOUBLE PRECISION A(NN), PIVOT(N), B(N), AFLAG(8) C INTEGER SNR(NN), RNR(NN1), HA(IHA,11), IFLAG(10) C INTEGER N, Z, NN, NN1, IHA, IFAIL C C These two versions can be used on many other computers C also. However some alterations may be needed and/or may C ensure greater efficiency of the performance of the C subroutine. For example, it will be much more efficient C to declare arrays SNR, RNR and (if possible) HA as C INTEGER*2 arrays on some IBM installations. C C C 3. Method C C Gaussian elimination is used in the decomposition of C matrix A. Pivotal interchanges are implemented as an C attempt to preserve both the stability of the C computations and the sparsity of the coefficient C matrix. Two triangular matrices L and U are computed so C that LU=PAQ (where P and Q are permutation matrices). C The right hand side vector b is also modified so that C vector c = L**(-1)*P*b is computed on successful exit. C In this way matrix L is sometimes not needed in the C further computations and may be removed if this is so. C C 4. Parameters of the subroutine C C C N - INTEGER. On entry N must contain the number of C equations in the system Ax=b. Unchanged on C exit. C C Z - INTEGER. On entry Z must contain the number of C non-zero elements in the coefficient matrix A C of the system Ax=b. Unchanged on exit. C C A - REAL (in the single precision version Y12MCE) C or DOUBLE PRECISION (in the double precision C version Y12MCF) array of length NN (see below). C On entry the first Z locations of array A must C contain the non-zero elements of the C coefficient matrix A of the system Ax = b C ordered by rows (in the preceding call of C Y12MB). The content of array A is modified by C subroutine Y12MC. On successful exit array A C will contain the non-zero elements of the upper C triangular matrix U (without the diagonal C elements of matrix U which can be found in C array PIVOT, see below) and sometimes also the C non-zero elements of the lower triangular C matrix L (without the diagonal elements which C are not stored). C C SNR - INTEGER array of length NN (see below). On C entry SNR(j), j = 1(1)Z, must contain the C column number of the non-zero element stored in C A(j). This is accomplished by subroutine C Y12MB. The content of array SNR is modified by C subroutine Y12MC. On successful exit array SNR C will contain the column numbers of the non-zero C elements of the upper triangular matrix U C (without the column numbers of the diagonal C elements of matrix U) and sometimes also the C column numbers of the non-zero elements of the C lower triangular matrix L (without the column C numbers of the diagonal elements of L). C C NN - INTEGER. On entry NN must contain the length of C array A and SNR. Restriction: NN.ge.2*Z. C Recommended value: 2*Z .le.NN.le. 3*Z if the C non-zero elements of matrix L will be removed C by subroutine Y12MC (i.e. if Y12MC is called C with IFLAG(5) = 1) or 3*Z.le. NN.le.5*Z if the C non-zero elements of matrix L will be kept by C subroutine Y12MC (i.e. if Y12MC is called with C IFLAG(5) = 2). Unchanged on exit. C C RNR - INTEGER array of length NN1 (see below). On C entry the first Z locations of array RNR must C contain the row numbers of the non-zero C elements of the coefficient matrix A of the C system Ax = b so that the row numbers of the C non-zero elements in the first column of matrix C A are located before the row numbers of the C non-zero elements in the second column of C matrix A, the row numbers of the non-zero C elements in the second column of matrix A are C located before the row numbers of the non-zero C elements in the third column of matrix A and so C on (this ordereing is prepared by subroutine C Y12MB). The content of array RNR is modified by C subroutine Y12MC. On successful exit all C components of array RNR will normally be zero. C C NN1 - INTEGER. On entry NN1 must contain the length C of the array RNR. Restriction: NN1.ge.Z. C Recommended value: 2*Z.le.NN1.le.3*Z. C Unchanged on exit. C C PIVOT - REAL (in the single precision version Y12MCE) C or DOUBLE PRECISION (in the double precision C version Y12MCF) array of length N. The content C of array PIVOT is modified by subroutine Y12MC. C On successful exit array PIVOT will contain the C pivotal elements ( the diagonal elements of C matrix U). This means that a small element (or C small elements) in array PIVOT on exit may C indicate numerical singularity of the C coefficient matrix A. Note that the smallest in C absolute value element in array PIVOT is stored C in AFLAG(8), see below. C C B - REAL (in the single precision version Y12MCE) C or DOUBLE PRECISION (in the double precision C version Y12MCF) array of length N. The right C hand side vector b of the system Ax = b must be C stored in B on entry. The content of array B is C modified by subroutine Y12MC. On successful C exit array B will contain the components of C vector c = L**(-1)*P*b. C C HA - INTEGER two-dimensional array. The length of C the first dimension of HA is IHA (see below). C The length of the second dimension of HA is 11. C The content of columns 1, 3, 4, 6, 7, 8 and 11 C is prepared by Y12MB (and thus they should not C be altered between calls of Y12MB and Y12MC). C The content of array HA is modified by C subroutine Y12MC. The non-zero elements C (without the diagonal elements) of row i in C matrix L (when this matrix is stored) are C between HA(i,1) and HA(i,2) -1. If the non-zero C elements of matrix L are not stored then C HA(i,1)=HA(i,2), i=1(1)N. The non-zero elements C of row i in matrix U (without the diagonal C elements) are stored between HA(i,2) and C HA(i,3). Information concerning the row and C column interchanges is stored in HA(i,7) and C HA(i,8), i=1(1)N-1. Information about the C largest number of non-zero elements found in C row i / column j during any stage of the C elimination is stored (when IFLAG(4)=1 only) C in HA(i,9)/HA(j,10). The other information C stored in array HA is not used in the further C computations. C C IHA - INTEGER. On entry IHA must contain the length C of the first dimension of array HA. C Restriction: IHA C .ge.N. Unchanged on exit. C C AFLAG - REAL (in the single precision version Y12MCE) C or DOUBLE PRECISION (in the double precision C version Y12MCF) array of length 8. The content C of the array can be described as follows: C C AFLAG(1) - Stability factor. An element can be C chosen as pivotal element only if C this element is larger (in absolute C value) than the absolute value of C the largest element in its row C divided by AFLAG(1). On entry C AFLAG(1) should contain a real C number larger than 1.0. If this is C not so then the subroutine sets C AFLAG(1) = 1.0005. Recommended C values: AFLAG(1) ranging from 4.0 C to 16.0. Unchanged on exit (when C correctly initialized). C C AFLAG(2) - Drop-tolerance. An element which in C the process of the computations C becomes smaller (in absolute value) C than the drop-tolerance is removed C from array A (and its row and C column numbers are removed from C arrays RNR and SNR). On entry C AFLAG(2) should contain some C positive small number or zero. C Recommended value AFLAG(2) = C 1.0E-12 when matrix A is not too C badly scaled and the magnitude of C the elements is 1, otherwise C smaller values of AFLAG(2) should C be used. Unchanged on exit. C C AFLAG(3) - The subroutine will stop the C computation when the growth factor C (parameter AFLAG(5), see below) C becomes larger than AFLAG(3). On C entry AFLAG(3) should contain a C large positive number. If C AFLAG(3)<1.0e5 then the subroutine C sets AFLAG(3)=1.0e5. Recommended C value AFLAG(3)=1.0E6. Unchanged C on exit (when correctly C initialized). C C AFLAG(4) - The subroutine will stop the C computation when the absolute value C of a current pivotal element is C smaller than AFLAG(4)*AFLAG(6) C (parameter AFLAG(6) is described C below). On entry AFLAG(4) must C contain a small non-negative C number. If AFLAG(4)<0 on entry C then the subroutine sets C AFLAG(4)=-AFLAG(4). Recommended C value AFLAG(4)=1.0E-12. Unchanged C on exit (when correctly C initialized). C C AFLAG(5) - Growth factor. The content of C parameter AFLAG(5) is modified by C subroutine Y12MC. After each stage C of the Gaussian elimination C subroutine Y12MC sets AFLAG(5) = C AFLAG(7)/AFLAG(6) (parameters C AFLAG(6) and AFLAG(7) are described C below). On exit large values of C parameters AFLAG(5) indicate that C an appreciable error in the C computed solution is possible. In C an extreme case , where C AFLAG(5)>AFLAG(3) the subroutine C will terminate the computations in C an attempt to prevent overflow (and C IFAIL will be set equal to 4). C C AFLAG(6) - On entry AFLAG(6) is equal to the C largest element in the coefficient C matrix A of the system Ax = b (set C by subroutine Y12MB). Unchanged on C exit. C C AFLAG(7) - On exit the largest (in absolute C value) element found during any C stage of the elimination will be C stored in AFLAG(7). The content of C parameter AFLAG(7) is modified by C subroutine Y12MC. C C AFLAG(8) - On entry the minimal (in absolute C value) pivotal element will be C stored in AFLAG(8). Small values of C AFLAG(8) indicate numerical C singularity of the coefficient C matrix A. We advise the user to C check this parameter on exit from C the calculation very carefully. The C content of parameter AFLAG(8) is C modified by the subroutine Y12MC. C C IFLAG - INTEGER array of length 10. The content of this C array is modified by subroutine Y12MC. The C contents of the components of this array can be C described as follows. C C IFLAG(1) - This parameter is used in connection C with the error diagnostics. The user C should set IFLAG(1).ge.0 before the C call of package Y12M (i.e. before C the first call of a subroutine of C this package). IFLAG(1) is used in C the error checks by the subroutines C and should not be altered by the C user between any two successive C calls of subroutines of the package. C IFLAG(1) will be equal to -2 on C successful exit from subroutine C Y12MC. Thus, the content of IFLAG(1) C is modified by subroutine Y12MC. C C IFLAG(2) - On entry IFLAG(2) must contain some C positive integer smaller than N. We C recommend IFLAG(2).le.3. If IFLAG(3) C = 0 then this parameter is ignored C by subroutine Y12MC. If C IFLAG((3).ge.0 then the pivotal C search at any stage of the C elimination (except possibly some of C the last stages) is carried out in C the IFLAG(2) rows which have least C numbers of non-zero elements. C Unchanged on exit. C C IFLAG(3) - On entry IFLAG(3) must contain 0, 1 C or 2. For general pivotal search C IFLAG(3) should be set equal to 1. C If IFLAG(3) = 2 then only diagonal C elements of the coefficient matrix A C will be selected as pivotal C elements. If IFLAG(3) = 0 then the C Gaussian elimination will be carried C out without any pivoting. C IFLAG(3)=0 or IFLAG(3)=2 (i.e. one C of the special pivotal strategies is C to be applied) should be used very C carefully because the error C diagnostics algorithm may not trace C all errors in this case. For C example, if the user attempts to use C IFLAG(3)=0 for matrix A which C contains zero elements on the main C diagonal, then the run will often be C stopped because a division by zero C occurs. Unchanged on exit. C C IFLAG(4) - On entry IFLAG(4) must contain 0, 1 C or 2. IFLAG(4) = 0 is the best C choice when (i) only one system is C to be solved, (ii) the first system C of a sequence of systems with the C same matrix (Ax = b1, Ax = b2, ..., C Ax = bp) is to be solved, (iii) when C any system in a sequence of systems C whose matrices are of different C structure is to be solved. IFLAG(4) C = 1 is the best choice when the C first system of a sequence of C systems whose matrices are of the C same structure is to be solved. In C this case IFLAG(4) = 2 is to be used C in the solution of any system after C the first one. Unchanged on exit. C C IFLAG(5) - On entry IFLAG(5) must contain 1 or C 2. If IFLAG(5) = 1 then the C non-zero elements of matrix L will C be removed. If IFLAG(5) = 2 then the C non-zero elements of matrix L will C be stored. Unchanged on exit. C C IFLAG(6) - On successful exit IFLAG(6) will be C equal to the number of "garbage" C collections in the row ordered list. C If IFLAG(6) is large then it is C better to choose a larger value of C NN with next calls of subroutine C Y12MC with the same or a similar C matrix A. This will lead to a C reduction in the computing time. C The content of IFLAG(6) is modified C by the subroutine Y12MC. C C IFLAG(7) - On successful exit IFLAG(7) will be C equal to the number of "garbage" C collections in the column ordered C list. If IFLAG(7) is large then it C is better to choose a larger value C of NN1 in the next calls of C subroutine Y12MC with the same or a C similar matrix A. This will lead to C a reduction in the computing time. C The content of IFLAG(7) is modified C by subroutine Y12MC. C C IFLAG(8) - On successful exit IFLAG(8) will be C equal to the maximal number of C non-zero elements kept in array A at C any stage of the Gaussian C elimination. If IFLAG(8) is much C smaller then NN (or NN1) then the C length NN (or NN1) can be chosen C smaller in next calls of subroutine C Y12MC with the same or a similar C matrix A. This will lead to a C reduction of the storage needed. The C content of IFLAG(8) is modified by C subroutine Y12MC. C C IFLAG(9) - The content of IFLAG(3) is modified C by subroutine Y12MC when IFLAG(4) = C 1 and ignored otherwise. The C minimal length NN1 such that Y12MC C can solve systems whose matrices are C of the same structure without C "garbage" collections in the column C ordered list and "movings" of C columns at the end of the column C ordered list is stored in IFLAG(9) C after the solution of the first C system in the sequence (with C IFLAG(4) = 1). C C IFLAG(10)- The content of IFLAG(10) is modified C by the subroutine Y12MC when C IFLAG(4) = 1 and ignored otherwise. C The minimal length NN such that C subroutine Y12MC can solve systems C whose matrices are of the same C structure without "garbage" C collections in the row ordered list C and "movings" of rows to the end of C the row ordered list is stored in C IFLAG(10) after the solution of the C first system in the sequence (with C IFLAG(4) = 1). C C IFAIL - Error diagnostics parameter. The content of C parameter IFAIL is modified by subroutine C Y12MC. On exit IFAIL = 0 if the subroutine has C not detected any error. Positive values of C IFAIL on exit show that some error has been C detected by the subroutine. Many of the error C diagnostics are common for all subroutines in C the package. Therefore the error diagnostics C are listed in a separate section, Section 7, of C this book. We advise the user to check the C value of this parameter on exit. C C 5. Error diagnostics C C Error diagnostics are given by positive values of the C parameter IFAIL (see above). We advise the user to C check carefully the value of this parameter on exit. C The error messages are listed in Section 7. C C 6. Auxiliary subroutines C C None. C C 7. Timing C C The time taken depends on the order of the matrix C (parameter N), the number of the non-zero elements in C the matrix (parameter Z), the magnitude of the non-zero C elements and their distribution in the matrix. C C 8. Storage C C There are no internally declared arrays. C C 9. Accuracy C C It is difficult to evaluate the accuracy of the C computed solution. Large values of parameter AFLAG(5), C the growth factor, indicate unstable computations C during the decomposition stage. Small values of C parameter AFLAG(8) can be considered as a signal for C numerical singularity. We must emphasize here that C normally much more reliable evaluations of the accuracy C achieved can be found by the use of iterative C refinement, i.e. by the use of subroutine Y12MF. By the C use of the latter subroutine the computing time will C often be reduced too. C C 10. Some remarks C C Remark 1 The values on entry of Y12MC are the same as C the values on exit of Y12MB for many C parameters. Therefore the content of N, Z, C A, SNR, NN, RNR, NN1, columns 1, 3, 4, 6, 7, C 8 and 11 of HA, AFLAG(6), IFLAG(1), IFLAG(4) C and IFAIL should not be changed between C calls of Y12MB and Y12MC. C C Remark 2 The call of Y12MC is normally followed by C one or several calls of Y12MD. The content C of N, A, SNR, NN, B, PIVOT, columns 1, 2, 3, C 7 and 8 of HA, IHA, AFLAG, IFLAG(1), C IFLAG(3), IFLAG(4) and IFAIL should not be C altered between calls of Y12MC and Y12MD. C C Remark 3 If IFAIL > 0 on exit from Y12MC there is no C justification for calling Y12MD. Therefore C the computations should be terminated and C the user should investigate (using the error C diagnostics given in Section 7) why Y12MC C has assigned a positive value to IFAIL. C ================================= C Documentation of subroutine Y12MD C ================================= C C 1. Purpose C C Y12MD solves sparse systems of linear algebraic C equations. C C 2. Calling sequence and declaration of the parameters C C The subroutine is written in FORTRAN and it has been C extensively tested with the FOR and FTN compilers on C the UNIVAC 1100/82 computer at the Regional Computing C Centre at the University of Copenhagen (RECKU). Many C examples have been run on an IBM 3033 computer at the C Northern Europe University Computing Centre (NEUCC) and C on a CDC Cyber 173 computer at the Regional Computing C Centre at the University of Aarhus (RECAU). Two C different versions are available: a single precision C version named Y12MDE and a double precision version C named Y12MDF. The calls of these two versions and the C declarations of the parameters are as follows. C C A). Single precision version: Y12MDE C C SUBROUTINE Y12MDE(N, A, NN, B, PIVOT, SNR, C 1 HA, IHA,, IFLAG, IFAIL) C REAL A(NN), PIVOT(N), B(N) C INTEGER SNR(NN), HA(IHA,11), IFLAG(10) C INTEGER N, NN, IHA, IFAIL C C B). Double precision version: Y12MDF C C SUBROUTINE Y12MDF(N, A, NN, B, PIVOT, SNR, C 1 HA, IHA,, IFLAG, IFAIL) C DOUBLE PRECISION A(NN), PIVOT(N), B(N) C INTEGER SNR(NN), HA(IHA,11), IFLAG(10) C INTEGER N, NN, IHA, IFAIL C C These two versions can be used on many others computers C also. However some alterations may be needed and/or may C ensure greater efficiency of the performance of the C subroutine. For example, it will be much more efficient C to declare arrays SNR and (if possible) HA as INTEGER*2 C arrays on some IBM installations. C C C 3. Method C C The LU decomposition (or only the upper triangular C matrix U if the vector c = L**(-1)*P*b is already C computed) is used to obtain an approximation to the C solution vector x of the system Ax = b. The LU C decomposition (or only matrix U ) must be available on C entry. This means that if a system with a new matrix is C solved then the call of subroutine Y12MD must be C preceded by a call of Y12MC and the contents of some C parameters and arrays should not be altered between the C calls of Y12MC and Y12MD. If a system with the same C matrix has already been solved, then Y12MD only should C be called (but the contents of some parameters and C arrays used in the preceding call of Y12MD should not C be altered). C C 4. Parameters of the subroutine C C N - INTEGER. On entry N must contain the number of C equations in the system Ax=b. Unchanged on C exit. C C A - REAL (in the single precision version Y12MDE) C or DOUBLE PRECISION (in the double precision C version Y12MDF) array of length NN (see below). C On entry array A must contain the non-zero C elements of the upper triangular matrix U C (without the diagonal elements) and sometimes C also the non-zero elements of the lower C triangular matrix L under the diagonal (these C elements are calculated by the subroutine C Y12MC, the non-zero elements of L are stored by C Y12MC only when this subroutine is called with C IFLAG(5) = 2). The content of the array A is C not modified by subroutine Y12MD. C C NN - INTEGER. On entry NN must contain the length of C the arrays A and SNR. Restriction: NN > Z. C Recommended value: 2*Z < NN < 3*Z if the C non-zero elements of matrix L will be removed C during the decomposition (i.e. if subroutine C Y12MC is called with IFLAG(5) = 1) and 3*Z < NN C < 5*Z if the non-zero elements of matrix L will C be stored during the decomposition (i.e. if the C subroutine Y12MC is called with IFLAG(5) = 2). C Unchanged on exit. C C B - REAL (in the single precision version Y12MDE) C or DOUBLE PRECISION (in the double precision C version Y12MDF) array of length N. If C subroutine Y12MC has been called in the C solution of system Ax = b (i.e. if C IFLAG(5)<3), then on entry array B must contain C vector c = L**(-1)*P*b (computed by Y12MC). If C subroutine Y12MC has not been called in the C solution of system Ax = b (i.e a system with C the same matrix has been solved before solving C Ax = b and the old LU decomposition is used by C setting IFLAG(5) = 3, see below) the array B C must contain the right hand side vector b on C entry. The content of the array B is modified C by the subroutine Y12MD. On successful exit an C approximation to the true solution x will be C found in B. C C PIVOT - REAL (in the single precision version Y12MDE) C or DOUBLE PRECISION (in the double precision C version Y12MDF) array of length N. On entry the C diagonal elements of the upper triangular C matrix U must be stored in array PIVOT (these C elements are calculated and stored in PIVOT by C subroutine Y12MC). Unchanged on exit. C C SNR - INTEGER array of length NN (see below). On C entry the column numbers of the non-zero C elements of the upper triangular matrix U C (without the column numbers of the non-zero C elements on the diagonal) and sometimes also C the column numbers of the non-zero elements of C the lower triangular matrix L (without the C column numbers of the non-zero elements on the C diagonal) must be stored in array SNR. This C structure is prepared by Y12MC. The content of C array SNR is not modified by the subroutine C Y12MD. C C HA - INTEGER two-dimesional array. The length of C the first dimension of HA is IHA (see below). C The length of the second dimension of HA is 11. C The content of columns 4, 5, 6, 9, 10 and 11 is C ignored by subroutine Y12MD. The content of C columns 1, 2, 3, 7 and 8 is stored by C subroutine Y12MC and should not be altered C between calls of Y12MC and Y12MD. Unchanged on C exit. C C IHA - INTEGER. On entry IHA must contain the length C of the first dimension of array HA. C Restriction: IHA .ge. N. Unchanged on exit. C C IFLAG - INTEGER array of length 10. The contents of all C locations in this array except IFLAG(1), C IFLAG(3), IFLAG(4) and IFLAG(5) are ignored by C subroutine Y12MD. The user should set IFLAG(1) C .ge. 0 before the call of package Y12M (i.e. C before the first call of a subroutine of this C package). IFLAG(1) is used in the error checks C by the subroutines and should not be altered by C the user between any two successive calls of C subroutines of the package. IFLAG(1) is equal C to -2 both on entry and on successful exit C from subroutine Y12MD. Thus, the content of C IFLAG(1) is not modified by subroutine Y12MD. C The same values of IFLAG(3) and IFLAG(4) as C those in the preceding call of Y12MC should be C used. If subroutine Y12MC has been called in C the solution of the system Ax = b, then on C entry IFLAG(5) must have the same value as on C entry of Y12MC. If subroutine Y12MC has not C been called in the solution of the system Ax = C b (i.e. a system with the same matrix has been C solved before solving Ax = b and the LU C decomposition of matrix A is thus available) C then on entry IFLAG(5) must contain the number C 3. The content of array IFLAG is unchanged on C exit. C C IFAIL - Error diagnostics parameter. The content of C parameter IFAIL is modified by subroutine C Y12MA. On exit IFAIL = 0 if the subroutine has C not detected any error. Positive values of C IFAIL on exit show that some error has been C detected by the subroutine. Many of the error C diagnostics are common for all subroutines in C the package. Therefore the error diagnostics C are listed in a separate section, Section 7, of C this book. We advise the user to check the C value of this parameter on exit. C C 5. Error diagnostics C C Error diagnostics are given by positive values of the C parameter IFAIL (see above). We advise the user to C check carefully the value of this parameter on exit. C The error messages are listed in Section 7. C C 6. Auxiliary subroutines C C None. C C 7. Timing C C The time taken depends on the order of the matrix C (parameter N) and the number of the non-zero elements C in the LU-decomposition of matrix A. C C 8. Storage C C There are no internally declared arrays. C C 9. Accuracy C C It is difficult to evaluate the accuracy of the C computed solution. Large values of parameter AFLAG(5), C the growth factor, indicate unstable computations C during the decomposition stage. Small values of C parameter AFLAG(8), the minimal pivotal element, can be C considered as a signal for numerical singularity. We C must emphasize here that normally much more reliable C evaluations of the accuracy achieved can be found by C the use of iterative refinement, i.e. by the use of C subroutine Y12MF. By the use of the latter subroutine C the computing time will often be reduced too. However, C the storage requirements may sometimes be increased. C C 10. Some remarks C C Remark 1 The content of N, A, SNR, NN, columns 1, 2, 3, C 7 and 8 of HA, IHA, IFLAG(3), IFLAG(4) and C IFAIL should not be altered between calls of C Y12MC and Y12MD. C C Remark 2 If the LU decomposition of matrix A is C available (i.e. a system with matrix A has C already been solved) then Y12MB and Y12MC C should not be called. The user should only C assign the new right hand side vector to array C B, set IFLAG(5) = 3 and call Y12MD. C ================================= C Documentation of subroutine Y12MF C ================================= C C 1. Purpose C C Large and sparse systems of linear algebraic equations C with real coefficients are solved by the use of C Gaussian elimination and sparse matrix technique. C Iterative refinement is applied in order to improve the C accuracy. C C 2. Calling sequence and declaration of the parameters C C The subroutine is written in FORTRAN and has been C extensively tested with the FOR and FTN compilers on C the UNIVAC 1100/82 computer at the Regional Computing C Centre at the University of Copenhagen (RECKU). Many C examples have been run on an IBM 3033 computer at the C Northern Europe University Computing Centre (NEUCC) and C on a CDC Cyber 173 computer at the Regional Computing C Centre at the University of Aarhus (RECAU). Only a C single precision version, Y12MFE, of this subroutine is C available. The call of the subroutine and the C declarations of the parameters are as follows. C C SUBROUTINE Y12MFE(N, A, SNR, NN, RNR, NN1, A1, C SN, NZ, C 1 HA, IHA, B, B1, X, Y, AFLAG, C IFLAG, IFAIL) C REAL A(NN), B(N), B1(N), X(N), Y(N), A1(NZ), C AFLAG(11) C INTEGER SNR(NN), RNR(NN1), HA(IHA,13), SN(NZ), C IFLAG(12) C INTEGER N, NN, NN1, NZ, IHA, IFAIL C C This subroutine can be used on many other computers C also. However, some alterations may be needed and/or C may ensure greater efficiency of the performance of the C subroutine. For example, it will be much more C efficient to declare arrays SNR, RNR, SN and (if C possible) HA as INTEGER*2 arrays on some IBM C instalations. Note too that the subroutine accumulates C some inner products in double precision and then rounds C them to single precision. Therefore if the use of C double precision is very expensive (in comparison with C the use of single precision) on the computer at the C user's disposal then it will not be efficient to use C subroutine Y12MF. The single precision versions of the C other subroutines should be used in this situation. C C 3. Method C C Consider the system Ax = b, where matrix A is sparse. C The non-zero elements of matrix A are ordered by rows C (subroutine Y12MB is called to perform this operation) C and then factorized into two triangular matrices, an C upper triangular matrix U and a unit lower triangular C matrix L (subroutine Y12MC is called to calculate the C non-zero elements of U amd L). The factorized system is C solved (subroutine Y12MD is called to calculate an C approximation x1 to the solution vector x) and an C attempt to improve the first solution by iterative C refinement is carried out in the following way: C C r(k-1) = b - A*x(k-1) C C d(k-1) = Q*U**(-1)*L(-1)*P*r(k-1) C C x(k) = X(k-1) + d(k-1) C C where k = 2, 3, ..., p. C C Normally, this process will improve the accuracy of the C first solution x1. If the process is not convergent or C the convegence is very slow, then information about the C trouble can be obtained by checking the parameters C AFLAG(10) (the max-norm of the last residual vector C r(p-1) computed by Y12MF) and AFLAG(9) (the max-norm of C the last correction vector, d(p-1), computed by Y12MF). C We advise the user to check carefully these parameters C on exit. Large values of these parameters show that the C computed solution is not very accurate. C C 4. Parameters of the subroutine C C N - INTEGER. On entry N must contain the number of C equations in the system Ax=b. Unchanged on C exit. C C A - REAL array of length NN (see below). If the C coefficient matrix must be factorized (in this C case IFLAG(5) = 2), then on entry the first NZ C locations of array A must contain the non-zero C elements of matrix A, the order of the elements C can be arbitrary. If the factorization of C matrix A is available (i.e. a system with the C same matrix has been solved in a previous call C of Y12MF) then on entry array A must contain C the non-zero elements of U and L (IFLAG(5) = 3 C should also be assigned in this case). The C content of array A is modified by Y12MF when C the coefficient matrix has to be factorized. C Unchanged on exit when an old LU factorization C is used. The content of array A should not be C altered between successive calls of Y12MF for C the solution of systems with the same matrices. C C SNR - INTEGER array of the length NN (see below). If C the coefficient matrix A must be factorized C then on entry the first NZ positions of array A C must contain the column numbers of the non-zero C elements of matrix A (ordered as the non-zero C elements in array A). The content of array SNR C is modified by the subroutine Y12MF in this C case. If the LU factorization of matrix A is C available then on entry array SNR must contain C the column numbers of the non-zero elements of C U and L. The content of array SNR is unchanged C on exit in this case. The content of array SNR C should not be altered between successive calls C of Y12MF for the solution of systems with the C same matrices. C C NN - INTEGER. On entry NN must contain the length of C array A and SNR. Restriction: NN .ge. 2*Z. C Recommended value: 2*Z .le. NN .le. 3*Z C Unchanged on exit. C C RNR - INTEGER array of length NN1 (see below). If the C coefficient matrix A must be factorized, then C on entry the first NZ positions of array RNR C must contain the row numbers of the non-zero C elements of matrix A (ordered as the non-zero C elements in array A). The content of array RNR C is modified by the subroutine in this case. If C the LU factorization of matrix A is available C then the content of array RNR is ignored by the C subroutine Y12MF. C C NN1 - INTEGER. On entry NN1 must contain the length C of the array RNR. Restriction: NN1 .ge. Z. C Recommended value: 1.5*Z .le. NN1 .le. 2*Z. C Unchanged on exit. C C A1 - REAL array of length NZ (see below). The C content of array A1 is modified by subroutine C Y12MF when the LU factorization of the C coefficient matrix A is not available. C Subroutine Y12MF copies the first NZ locations C of array A into array A1 in this case (after C the internal call of Y12MB; this means that C array A1 contains the non-zero elements of the C coefficient matrix A ordered by rows on exit). C If the LU factorization of the coefficient C matrix A is available then the content of array C A1 is unchanged on exit. The content of array C A1 should not be altered between successve C calls of Y12MF in the solution of systems with C the same matrices. C C SN - INTEGER array of length NZ (see below). The C content of array SN is modified by subroutine C Y12MF when the LU factorization of the C coefficient matrix A is not available. C Subroutine Y12MF copies the first NZ locations C of array SNR into array SN in this case (after C the internal call of Y12MB; this means that C array SN contains the column numbers of the C non-zero elements ordered by rows on exit). If C the LU factorization is available then the C content of array SN is unchanged on exit. The C content of array SN should not be altered C between successive calls of Y12MF in the C solution of systems with the same matrices. C C NZ - INTEGER. On entry NZ must contain the number of C non-zero elements in the coefficient matrix A C of the system Ax = b. Unchanged on exit. C C HA - INTEGER two-dimensional array. The length of C the first dimension of HA is IHA (see below). C The length of the second dimension of HA is 13. C If a new decomposition should be calculated, C then subroutine Y12MF modifies the contents of C HA. The contents on exit of some of the columns C of array HA are described in the documentation C of subroutine Y12MC. The last two columns C (12'th and 13'th) contain on exit information C about the row starts and the row ends in the C original matrix (after the non-zero elements of C this matrix have been ordered by rows); this C means that the non-zero elements in row i of C the coefficient matrix A are located between C positions HA(i,12) and HA(i,13) in array A1 C (the column numbers of the non-zero elements in C row i are also located between positions C HA(i,12) and HA(i,13) in array SN). If the C matrix of the system has the same structure as C the matrix of a system which is already solved, C then subroutine Y12MF will use the information C stored in columns 7, 8, 9 and 10 of HA during C the previous call (therefore this information C should not be altered). If the LU decomposition C of the coefficient matrix is available, then C the contents of array HA are unchanged on exit. C The contents of columns 1, 2, 3, 7, 8, 12 and C 13 of array HA should not be altered between C successive calls of Y12MF in the solution of C systems with the same matrix. C C IHA - INTEGER. On entry IHA must contain the length C of the first dimension of array HA. C Restriction: IHA .ge. N. Unchanged on exit. C C B - REAL array of length N. On entry the right-hand C side vector b of the system Ax=b must be C stored in array B. The content of array B is C modified by subroutine Y12MF. On successful C exit the components of the last correction C vector d(p-1) are stored in array B. C C B1 - REAL array of length N. The content of array B1 C is modified by subroutine Y12MF. The right hand C side vector of the system Ax = b is stored in C array B1 on successful exit. C C X - REAL array of length N. The content of array X C is modified by the subroutine Y12MF. On C successful exit the corrected solution vector C is stored in array X. C C Y - REAL array of length N. The content of array Y C is modified by subroutine Y12MF. On successful C exit array Y will contain the pivotal elements C (the diagonal elements of matrix U). This means C that a small element (or small elements) in C array Y on exit may indicate numerical C singularity of the coefficient matrix A. Note C that the smallest in absolute value element in C array Y is also stored in AFLAG(8), see below. C C AFLAG - REAL array of length 11. The content of the C array can be described as follows: C C AFLAG(1) - Stability factor. An element can be C chosen as pivotal element only if C this element is larger (in absolute C value) than the absolute value of C the largest element in its row C divided by AFLAG(1). On entry C AFLAG(1) should contain a real C number larger than 1.0. If this is C not so then the subroutine sets C AFLAG(1) = 1.0005. Recommended C values of AFLAG(1) ranging from 4.0 C to 16.0. Unchanged on exit (when C correctly initialized). C C AFLAG(2) - Drop-tolerance. An element which in C the process of the computations C becomes smaller(in absolute value) C than the drop-tolerance is removed C from array A (and its row and C column numbers are removed from C arrays RNR and SNR). Recommended C value: on entry AFLAG(2) should C be in the interval C (a*1.0E-4,a*1.0E-1) where a is the C magnitude of the elements. If the C magnitude a of the non-zero C elements is not known, then the C user should intialize the required C drop-tolerance multiplied by -1. C Let a(i) be the maximal element (in C absolute value) in row i ( i=1(1)N C ). Let a be the minimal number C among all a(i). If a negative value C of AFLAG(2) is assigned, the C subroutine computes a and uses a C drop-tolerance equal to C -AFLAG(2)*a. In the above C considerations it is assumed that C the coefficient matrix A is not too C badly scaled. If this is not so C then perform row scaling of the C coefficient matrix so that the C largest elements in the rows are of C magnitude 1 and choose AFLAG(2) in C the interval (1.0E-5,1.0E-3). C Unchanged on exit. C C AFLAG(3) - The subroutine will stop the C computation when the growth factor C (parameter AFLAG(5), see below) C becomes larger than AFLAG(3). On C entry AFLAG(3) should contain a C large positive number. If C AFLAG(3)<1.0E5 then the subroutine C sets AFLAG(3)=1.0E5. Recommended C value AFLAG(3)=1.0E16. Unchanged C on exit (when correctly C initialized). C C AFLAG(4) - The subroutine will stop the C computation when the absolute value C of a current pivotal element is C smaller than AFLAG(4)*AFLAG(6) C (parameter AFLAG(6) is described C below). On entry AFLAG(4) must C contain a small non-negative C number. If AFLAG(4)<0 then the C subroutine sets C AFLAG(4)=-AFLAG(4). Recommended C value AFLAG(4)=1.0E-12. Unchanged C on exit (when correctly C intialized). C C AFLAG(5) - Growth factor. The content of C parameter AFLAG(5) is modified by C subroutine Y12MF. After each stage C of the Gaussian elimination C subroutine Y12MF sets AFLAG(5) = C AFLAG(7)/AFLAG(6) (parameters C AFLAG(6) and AFLAG(7) are described C below). On exit large values of C parameters AFLAG(5) indicate that C an appreciable error in the C computed solution is possible. In C an extreme case, where C AFLAG(5)>AFLAG(3), the subroutine C will terminate the computations in C an attempt to prevent overflow. C C AFLAG(6) - On exit AFLAG(6) is equal to the C largest element in the coefficient C matrix A of the system Ax=b (set C by subroutine Y12MB). Unchanged on C exit. C C AFLAG(7) - On exit the largest (in absolute C value) element found during any C stage of the elimination will be C stored in AFLAG(7). The content of C parameter AFLAG(7) is modified by C subroutine Y12MF. C C AFLAG(8) - On exit the minimal (in absolute C value) pivotal element will be C stored in AFLAG(8). Small values of C AFLAG(8) indicate numerical C singularity of the coefficient C matrix A. We advise the user to C check this parameter on exit from C the calculation very carefully. The C content of parameter AFLAG(8) is C modified by the subroutine Y12MF. C C AFLAG(9) - The content of AFLAG(9) is modified C by the subroutine Y12MF. The C max-norm of the last correction C vector d(p-1) will be stored in C AFLAG(9) on successful exit. We C advise the user to check carefully C this parameter. C C AFLAG(10) - The content of AFLAG(10) is C modified by the subroutine Y12MF. C The max-norm of the last residual C vector r(p-1) will be stored in C AFLAG(10) on successful exit. We C advise the user to check carefully C this parameter. C C AFLAG(11) - The content of AFLAG(11) is C modified by the subroutine Y12MF. C On exit AFLAG(11) will contain the C max-norm of the corrected solution C vector. If the value of AFLAG(11) C is not to close to zero then C AFLAG(9)/AFLAG(11) will normally C give an excellent evaluation of the C relative error in the solution C vector. C C IFLAG - INTEGER array of length 12. The content of this C array can be described as follows: C C IFLAG(1) - This parameter is used as a work C space by subroutine Y12MF. C C IFLAG(2) - On entry IFLAG(2) must contain some C positive integer smaller than N. We C recommend IFLAG(2) .le. 3. If C IFLAG(3) = 0 then this parameter is C ignored by subroutine Y12MF. If C IFLAG((2) .ge. 0 then the pivotal C search at any stage of the C elimination (except possibly some of C the last stages) is carried out in C the IFLAG(2) rows which have least C number of non-zero elements. C Unchanged on exit. C C IFLAG(3) - On entry IFLAG(3) must contain 0, 1 C or 2. For general pivotal search C IFLAG(3) should be set equal to 1. C If IFLAG(3) = 2 then only diagonal C elements of the coefficient matrix A C can be selected as pivotal elements. C If IFLAG(3) = 0 then the Gaussian C elimination will be carried out C without any pivoting. IFLAG(3)=0 or C IFLAG(3)=2 (i.e. one of the special C pivotal strategies is to be applied) C should be used very carefully C because the error diagnostics C algorithm may not trace all errors C in this case. Unchanged on exit. C C IFLAG(4) - On entry IFLAG(4) must contain 0, 1 C or 2. IFLAG(4) = 0 is the best C choice when (i) only one system is C to be solved, (ii) the first system C of a sequence of systems with the C same matrix (Ax = b1, Ax = C b2,......Ax = bp) is to be solved, C (iii) when any system in a sequence C of systems whose matrices are of C different structure is to be solved. C IFLAG(4) = 1 is the best choice when C the first system of a sequence of C systems whose matrices are of the C same structure is to be solved. In C this case IFLAG(4) = 2 can be used C in the solution of any system after C the first one. Unchanged on exit. C C IFLAG(5) - If the LU factorization of the C coefficient matrix is not available, C then IFLAG(5) must be set to 2 on C entry. If the LU factorization of C the coefficient matrix is available, C then IFLAG(5) must be set to 3 on C entry. Unchanged on exit. C C IFLAG(6) - On successful exit IFLAG(6) will be C equal to the number of "garbage" C collections in the row ordered list. C If IFLAG(6) is large then it is C better to choose a larger value of C NN with next calls of subroutine C Y12MF with the same or similar C matrix A. This will lead to a C reduction in the computing time. C The content of IFLAG(6) is modified C by the subroutine Y12MF. C C IFLAG(7) - On successful exit IFLAG(7) will be C equal to the number of "garbage" C collections in the column ordered C list. If IFLAG(7) is large then it C is better to choose a larger value C of NN1 in the next calls of C subroutine Y12MF with the same or C similar matrix A. This will lead to C a reduction of the computing time. C The content of IFLAG(7) is modified C by subroutine Y12MF. C C IFLAG(8) - On successful exit IFLAG(8) will be C equal to the maximal number of C non-zero elements kept in array A at C any stage of the Gaussian C elimination. If IFLAG(8) is much C smaller than NN (or NN1) then the C length NN (or NN1) can be chosen C smaller in next calls of subroutine C Y12MF with the same or similar C matrix A. This will lead to a C reduction of the storage needed. The C content of IFLAG(8) is modified by C subroutine Y12MF. C C IFLAG(9) - The minimal length NN1 such that C Y12MF can solve systems whose C matrices are of the same structure C without "garbage" collections in the C column ordered list and "movings" of C columns at the end of the column C ordered list is stored in IFLAG(9) C after the solution of the first C system in the sequence (with C IFLAG(4)=1). The content of C IFLAG(9) is modified by subroutine C Y12MF when IFLAG(4) = 1 and ignored C otherwise. C C IFLAG(10) The minimal length NN such that C subroutine Y12MF can solve systems C whose matrices are of the same C structure without "garbage" C collections in the row ordered list C and "movings" of rows to the end of C the row ordered list is stored in C IFLAG(10) after the solution of the C first system in the sequence (with C IFLAG(4) = 1). The content of C IFLAG(10) is modified by the C subroutine Y12MF when IFLAG(4) = 1 C and ignored otherwise. C C IFLAG(11) The maximum allowed number of C iterations must be contained in C IFLAG(11) on entry. Restriction: C IFLAG(11)>1. Recommended value : C IFLAG(11)< 33. Unchanged on exit. C C IFLAG(12) The content of IFLAG(12) is modified C by the subroutine Y12MF. On exit C the number of iterations actually C performed will be stored in C IFLAG(12). C C IFAIL - Error diagnostic parameter. The content of C parameter IFAIL is modified by subroutine C Y12MF. On exit IFAIL = 0 if the subroutine has C not detected any error. Positive values of C IFAIL on exit show that some error has been C detected by the subroutine. Many of the error C diagnostics are common for all subroutines in C the package. Therefore the error diagnostics C are listed in a separate section, Section 7, of C this book. We advise the user to check the C value of this parameter on exit. C C 5. Error diagnostics C C Error diagnostics are given by positive values of C parameter IFAIL (see above). We advise the user to C check carefully the value of this parameter on exit. C The error messages are listed in Section 7. C C 6. Auxiliary subroutines C C Y12MF calls three other subroutines: Y12MB, Y12MC and C Y12MD. C C 7. Timing C C The time taken depends on the order of the matrix C (parameter N), the number of the non-zero elements in C the matrix (parameter Z), the magnitude of the non-zero C elements and their distribution in the matrix. C C 8. Storage C C There are no internally declared arrays. C C 9. Accuracy C C Normally the accuracy achieved can be estimated very C well by means of AFLAG(9) and AFLAG(10). Further C information about the accuracy can sometimes be C obtained by inspection of the growth factor AFLAG(5), C and the smallest pivotal element (in absolute value) C stored in AFLAG(8). C C 10. Some remarks C C Remark 1 Following the recommendations given in Section C 1.3.6 the user can write a subroutine (similar C to the subroutine Y12MA) where the recommended C values of the parameters AFLAG(1) to AFLAG(4), C IFLAG(2) to IFLAG(5) and IFLAG(11) are C initialized. If this is done then only N, NN, C NN1, IHA, A, SNR, RNR and B should be assigned C before the call of this new subroutine. C C Remark 2 The information stored in some of the arrays C should not be altered between successive calls C of subroutine Y12MF (see more details in C Section 1.3.6). C ================= C Error diagnostics C ================= C C IFAIL is the error diagnostics parameter. On exit from C each subroutine IFAIL = 0 if no error has been C detected. IFAIL > 0 indicates that some error has been C detected and the computations have been terminated C immediately after the detection of the error. The C errors corresponding to the different positive values C of IFAIL are listed below. C C IFAIL = 1 The coefficient matrix A is not C factorized, i.e. the call of subroutine C Y12MD was not preceded by a call of Y12MC C during the solution of Ax=b or during C the solution of the first system in a C sequence ( Ax1 = b1 , Ax2 = b2,.....,Axp = C bp) of systems with the same coefficient C matrix. This will work in all cases only C if the user sets IFLAG(1) .ge. 0 before C the call of package Y12M (i.e. before the C first call of a subroutine of this C package). C C IFAIL = 2 The coefficient matrix A is not ordered, C i.e. the call of subroutine Y12MC was not C preceded by a call of Y12MB. This will C work in all cases only if the user sets C IFLAG(1) .ge. 0 before the call of package C Y12M (i.e. before the first call of a C subroutine of this package). C C IFAIL = 3 A pivotal element abs(a(i,i;j)) < AFLAG(4) C * AFLAG(6) is selected. When AFLAG(4) is C sufficiently small this is an indication C that the coefficient matrix is numerically C singular. C C IFAIL = 4 AFLAG(5), the growth factor, is larger C than AFLAG(3). When AFLAG(3) is C sufficiently large this indicates that the C elements of the coefficient matrix A grow C so quickly during the factorization that C the continuation of the computation is not C justified. The choice of a smaller C stability factor, AFLAG(1), may give C better results in this case. C C IFAIL = 5 The length NN of arrays A and SNR is not C sufficient. Larger values of NN (and C possibly of NN1) should be used. C C IFAIL = 6 The length NN1 of array RNR is not C sufficient. Larger values of NN1 (and C possibly of NN) should be used. C C IFAIL = 7 A row without non-zero elements in its C active part is found during the C decomposition. If the drop-tolerance, C AFLAG(2), is sufficiently small, then C IFAIL = 7 indicates that the matrix is C numerically singular. If a large value of C the drop-tolerance AFLAG(2) is used and if C IFAIL = 7 on exit, this is not certain. A C run with a smaller value of AFLAG(2) C and/or a careful check of the parameters C AFLAG(8) and AFLAG(5) is recommended in C the latter case. C C IFAIL = 8 A column without non-zero elements in its C active part is found during the C decomposition. If the drop-tolerance, C AFLAG(2), is sufficiently small, then C IFAIL = 8 indicates that the matrix is C numerically singular. If a large value of C the drop-tolerance AFLAG(2) is used and if C IFAIL = 8 on exit, this is not certain. A C run with a smaller value of AFLAG(2) C and/or a careful check of the parameters C AFLAG(8) and AFLAG(5) is recommended in C the latter case. C C IFAIL = 9 A pivotal element is missing. This may C occur if AFLAG(2) > 0 and IFLAG(4) = 2 C (i.e. some system after the first one in a C sequence of systems with the same C structure is solved using a positive value C for the drop-tolerance). The value of the C drop-tolerance AFLAG(2), should be C decreased and the coefficient matrix of C the system refactorized. This error may C also occur when one of the special pivotal C strategies (IFLAG(3)=0 or IFLAG(3)=2) is C used and the matrix is not suitable for C such a strategy. C C IFAIL = 10 Subroutine Y12MF is called with IFLAG(5) = C 1 (i.e. with a request to remove the C non-zero elements of the lower triangular C matrix L). IFLAG(5)=2 must be C initialized instead of IFLAG(5)=1. C C IFAIL = 11 The coefficient matrix A contains at least C two elements in the same position (i,j). C The input data should be examined C carefully in this case. C C IFAIL = 12 The number of equations in the system Ax=b C is smaller than 2 (i.e. N<2). The value C of N should be checked. C C IFAIL = 13 The number of non-zero elements of the C coefficient matrix is non-positive (i.e. C Z.le.0 ). The value of the parameter Z C (renamed NZ in Y12MF) should be checked. C C IFAIL = 14 The number of non-zero elements in the C coefficient matrix is smaller than the C number of equations (i.e. Z < N ). If C there is no mistake (i.e. if parameter Z, C renamed NZ in Y12MF, is correctly assigned C on entry) then the coefficient matrix is C structurally singular in this case. C C IFAIL = 15 The length IHA of the first dimension of C array HA is smaller than N. IHA.ge.N C should be assigned. C C IFAIL = 16 The value of parameter IFLAG(4) is not C assigned correctly. IFLAG(4) should be C equal to 0, 1 or 2. See more details in C the description of this parameter. C C IFAIL = 17 A row without non-zero elements has been C found in the coefficient matrix A of the C system before the Gaussian elimination is C initiated. Matrix A is structurally C singular. C C IFAIL = 18 A column without non-zero elements has C been found in the coefficient matrix A of C the system before the Gaussian elimination C is initiated. Matrix A is structurally C singular. C C IFAIL = 19 Parameter IFLAG(2) is smaller than 1. The C value of IFLAG(2) should be a positive C integer (IFLAG(2) = 3 is recommended). C C IFAIL = 20 Parameter IFLAG(3) is out of range. C IFLAG(3) should be equal to 0, 1 or 2. C C IFAIL = 21 Parameter IFLAG(5) is out of range. C IFLAG(5) should be equal to 1, 2 or 3 (but C when IFLAG(5) = 3 Y12MB and Y12MC should C not be called; see also the message for C IFAIL = 22 below). C C IFAIL = 22 Either subroutine Y12MB or subroutine C Y12MC is called with IFLAG(5) = 3. Each of C these subroutines should be called with C IFLAG(5) equal to 1 or 2. C C IFAIL = 23 The number of allowed iterations C (parameter IFLAG(11) when Y12MF is used) C is smaller than 2. IFLAG(11) .ge. 2 C should be assigned. C C IFAIL = 24 At least one element whose column number C is either larger than N or smaller than 1 C is found. C C IFAIL = 25 At least one element whose row number is C either larger than N or smaller than 1 is C found. C ========== C References C ========== C 1. Bjorck, A. - C "Methods for Sparse Linear Least-Squares C Problems". C In: "Sparse Matrix Computations" C (J.Bunch and D.Rose, eds.), pp.177-199. C Academic Press, New York, 1976. C C 2. Clasen, R.J. - C "Techniques for Automatic Tolerance in C Linear Programming", C Comm. ACM 9, pp. 802-803, 1966. C C 3. Cline, A.K., Moler, C.B., Stewart, G.W. and C Wilkinson, J.H. - C "An estimate for the condition number of a C matrix", C SIAM J. Numer. Anal. 16, 368-375, 1979. C C 4. Dongarra, J.J., Bunch, J.R., Moler, C.B. and C Stewart, G.W. - C "LINPACK User's Guide", C SIAM, Philadelphia, 1979. C C 5. Duff, I.S. - C "MA28 - a Set of FORTRAN Subroutines C for Sparse Unsymmetric Matrices", C Report No. R8730, A.E.R.E.,Harwell, C England, 1977. C C 6. Duff, I.S. and Reid, J.K. - C "Some Design Features of a Sparse Matrix C Code", C ACM Trans. Math. Software 5, pp. 18-35, C 1979. C C 7. Forsythe, G.E., Malcolm, M.A., and Moler, C.B. - C "Computer Methods for Mathematical C Computations", C Prentice-Hall, Englewood Cliffs, N.J., C 1977. C C 8. Forsythe, G.E. and Moler, C.B. - C "Computer Solution of Linear Algebraic C Equations", C Prentice-Hall, Englewood Cliffs, N.J., C 1967. C C 9. Gustavson, F.G. - C "Some Basic Techniques for Solving Sparse C Systems of Linear Equations". C In: "Sparse Matrices and Their C Applications", C (D.J. Rose and R.A. Willoughby, eds.), pp C 41-52, C Plenum Press, New York, 1972. C C 10. Gustavson, F.G. - C "Two Fast Algorithms for Sparse Matrices: C Multiplication and Permuted C Transposition", C ACM Trans. Math. Software, 4, pp. 250-269, C 1978. C C 11. Houbak, N. and Thomsen, P.G. - C "SPARKS - a FORTRAN Subroutine for C Solution C of Large Systems of Stiff ODE's with C Sparse Jacobians", C Report 79-02, Institute for Numerical C Analysis, C Technical University of Denmark, Lyngby, C Denmark, 1979. C C 12. Moler, C.B. - C "Three Research Problems in Numerical C Linear Algebra". C In: "Proceedings of the Symposia in C Applied Mathematics" C (G.H. Golub and J. Oliger, eds.), Vol. 22, C pp. 1-18, C American Mathematical Society, C Providence, Rhode Island, 1978. C C 13. NAG Library - C Fortran Manual, Mark 7, Vol. 3, 4, C Numerical Algorithms Group, C 7 Banbury Road, Oxford OX2 6NN, United C Kingdom. C C 14. Reid, J.K. - C "A Note on the Stability of Gaussian C Elimination", C J. Inst. Math. Appl., 8, pp. 374-375, C 1971. C C 15. Reid, J.K. - C "Fortran Subroutines for Handling Sparse C Linear Programming Bases", C Report R8269, A.E.R.E., Harwell, England, C 1976. C C 16. Schaumburg, K. and Wasniewski, J. - C "Use of a Semiexplicit Runge-Kutta C Integration Algorithm in a Spectroscopic C Problem" C Computers and Chemistry 2, pp. 19-25, C 1978. C C 17. Schaumburg, K., Wasniewski, J. and Zlatev, Z. - C "Solution of Ordinary Differential C Equations. Development of a Semiexplicit C Runge-Kutta Algorithm and Application to a C Spectroscopic Problem", C Computers and Chemistry, 3, pp. 57-63, C 1979. C C 18. Schaumburg, K., Wasniewski, J. and Zlatev, Z. - C "The Use of Sparse Matrix Technique in the C Numerical Integration of Stiff Systems of C Linear Ordinary Differential Equations", C Computers and Chemistry, 4, pp. 1-12, C 1980. C C 19. Stewart, G.W. - C "Introduction to Matrix Computations", C Academic Press, New York, 1973. C C 20. Tewarson, R.P. - C "Sparse Matrices", C Academic Press, New York, 1973. C C 21. Thomsen, P.G. - C "Numerical Solution of Large Systems with C Sparse Jacobians". C In: "Working Papers for the 1979 SIGNUM C Meeting on Numerical Ordinary Differential C Equations" (R.D. Skeel, ed.), C Computer Science Department, C University of Illinois at Urbana - C Champaign, Urbana, Illinois, 1979. C C 22. Wasniewski, J., Zlatev, Z. and Schaumburg, K. - C "A Multibanking Option of an Iterative C Refinement Subroutine". C In: "Conference Proceedings and Technical C Papers", C Spring Conference of Univac Users C Assotiation/Europe, Geneva, 1981. C C 23. Wilkinson, J.H. - C "Rounding Errors in Algebraic Processes", C Prentice-Hall, Englewood Cliffs, N.J., C 1963. C C 24. Wilkinson, J.H. - C "The Algebraic Eigenvalue Problem", C Oxford University Press, London, 1965. C C 25. Wilkinson, J.H and Reinsch, C. - C "Handbook for Automatic Computation", C Vol II, Linear Algebra, pp. 50-56, C Springer, Heidelberg, 1971. C C 26. Wolfe, P. - C "Error in the Solution of Linear C Programming Problems", C In: "Error in Digital Computation" C (L.B.Rall, ed.), Vol 2, pp. 271-284, C Wiley, New York, 1965. C C 27. Zlatev, Z. - C "Use of Iterative Refinement in the C Solution of Sparse Linear Systems", C Report 1/79, Institute of Mathematics and C Statistics, The Royal Veterinary and C Agricultural University, C Copenhagen, Denmark, 1979 C (to appear in SIAM J. Numer. Anal.). C C 28. Zlatev, Z. - C "On Some Pivotal Strategies in Gaussian C Elimination by Sparse Technique", C SIAM J. Numer. Anal. 17, pp. 18-30, 1980. C C 29. Zlatev, Z. - C "On Solving Some Large Linear Problems by C Direct Methods", C Report 111, Department of Computer C Science, University of Aarhus, Aarhus, C Denmark, 1980. C C 30. Zlatev, Z. - C "Modified Diagonally Implicit Runge-Kutta C Methods", C Report No. 112, Department of Computer C Science, University of Aarhus, Aarhus, C Denmark, 1980 C (to appear in SIAM Journal on Scientific C and Statistical Computing). C C 31. Zlatev, Z. - C "Comparison of Two Pivotal Strategies in C Sparse Plane Rotations", C Report 122, Department of Computer C Science, University of Aarhus, Aarhus, C Denmark, 1980. C (to appear in Computers and Mathematics C with Applications). C C 32. Zlatev, Z., Barker, V.A. and Thomsen, P.G. - C "SSLEST - a FORTRAN IV Subroutine for C Solving Sparse Systems of Linear Equations C (USER's GUIDE)", C Report 78-01, Institute for Numerical C Analysis, C Technical University of Denmark, Lyngby, C Denmark, 1978. C C 33. Zlatev, Z. and Nielsen, H.B. - C "Least - Squares Solution of Large Linear C Problems". C In: "Symposium i Anvendt Statistik 1980" C (A. Hoskuldsson, K. Conradsen, B. Sloth C Jensen and K.Esbensen, eds.), pp. 17-52. C NEUCC, Technical University of Denmark, C Lyngby, Denmark, 1980. C C 34. Zlatev, Z., Schaumburg, K. and Wasniewski, J. - C "Implementation of an Iterative Refinement C Option in a Code for Large and Sparse C Systems". C Computers and Chemistry, 4, pp. 87-99, C 1980. C C 35. Zlatev, Z. and Thomsen., P.G. - C "ST - a FORTRAN IV Subroutine for the C Solution of Large Systems of Linear C Algebraic Equations with Real Coefficients C by Use of Sparse Technique", C Report 76-05, Institute for Numerical C Analysis, C Technical University of Denmark, Lyngby, C Denmark, 1976. C C 36. Zlatev, Z. and Thomsen, P.G. - C "An Algorithm for the Solution of C Parabolic Partial Differential Equations C Based on Finite Element Discretization", C Report 77-09, Institute for Numerical C Analysis, C Technical University of Denmark, Lyngby, C Denmark, 1977. C C 37. Zlatev, Z. and Thomsen, P.G. - C "Application of Backward Differentiation C Methods to the Finite Element Solution of C Time Dependent Problems", C International Journal for Numerical C Methods in Engineering, 14, pp. 1051 - C 1061, 1979. C C 38. Zlatev, Z., Wasniewski, J. and Schaumburg, K. - C "Comparison of Two Algorithms for Solving C Large Linear Systems". C Report No 80/9, C Regional Computing Centre at the C University of Copenhagen, Vermundsgade 5, C DK-2100 Copenhagen, Denmark, 1980. C C 39. Zlatev, Z., Wasniewski, J. and Schaumburg, K. - C "Classification of the Systems of Ordinary C Differential Equations and Practical C Aspects in the Numerical Integration of C Large Systems", C Computers and Chemistry, 4, pp. 13-18, C 1980. C C 40. Zlatev, Z., Wasniewski, J., Schaumburg, K. - C "A Testing Scheme For Subroutines Solving C Large Linear Problems", C Report No 81/1, C Regional Computing Centre at the C University of Copenhagen, Vermundsgade 5, C DK-2100 Copenhagen, Denmark, 1981 C (to appear in Computers and Chemistry, 5, C 1981). Acknowledge-To: <NEUJW@vm.uni-c.dk>