2
SUBROUTINE SSDCGS (N, B, X, NELT, IA, JA, A, ISYM, ITOL, TOL,
3
+ ITMAX, ITER, ERR, IERR, IUNIT, RWORK, LENW, IWORK, LENIW)
4
C***BEGIN PROLOGUE SSDCGS
5
C***PURPOSE Diagonally Scaled CGS Sparse Ax=b Solver.
6
C Routine to solve a linear system Ax = b using the
7
C BiConjugate Gradient Squared method with diagonal scaling.
8
C***LIBRARY SLATEC (SLAP)
9
C***CATEGORY D2A4, D2B4
10
C***TYPE SINGLE PRECISION (SSDCGS-S, DSDCGS-D)
11
C***KEYWORDS ITERATIVE PRECONDITION, NON-SYMMETRIC LINEAR SYSTEM, SLAP,
13
C***AUTHOR Greenbaum, Anne, (Courant Institute)
14
C Seager, Mark K., (LLNL)
15
C Lawrence Livermore National Laboratory
17
C Livermore, CA 94550 (510) 423-3141
22
C INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX
23
C INTEGER ITER, IERR, IUNIT, LENW, IWORK(10), LENIW
24
C REAL B(N), X(N), A(NELT), TOL, ERR, RWORK(8*N)
26
C CALL SSDCGS(N, B, X, NELT, IA, JA, A, ISYM, ITOL, TOL,
27
C $ ITMAX, ITER, ERR, IERR, IUNIT, RWORK, LENW, IWORK, LENIW )
31
C Order of the Matrix.
33
C Right-hand side vector.
35
C On input X is your initial guess for solution vector.
36
C On output X is the final approximate solution.
38
C Number of Non-Zeros stored in A.
39
C IA :INOUT Integer IA(NELT).
40
C JA :INOUT Integer JA(NELT).
41
C A :INOUT Real A(NELT).
42
C These arrays should hold the matrix A in either the SLAP
43
C Triad format or the SLAP Column format. See "Description",
44
C below. If the SLAP Triad format is chosen it is changed
45
C internally to the SLAP Column format.
47
C Flag to indicate symmetric storage format.
48
C If ISYM=0, all non-zero entries of the matrix are stored.
49
C If ISYM=1, the matrix is symmetric, and only the upper
50
C or lower triangle of the matrix is stored.
52
C Flag to indicate type of convergence criterion.
53
C If ITOL=1, iteration stops when the 2-norm of the residual
54
C divided by the 2-norm of the right-hand side is less than TOL.
55
C This routine must calculate the residual from R = A*X - B.
56
C This is unnatural and hence expensive for this type of iter-
57
C ative method. ITOL=2 is *STRONGLY* recommended.
58
C If ITOL=2, iteration stops when the 2-norm of M-inv times the
59
C residual divided by the 2-norm of M-inv times the right hand
60
C side is less than TOL, where M-inv time a vector is the pre-
61
C conditioning step. This is the *NATURAL* stopping for this
62
C iterative method and is *STRONGLY* recommended.
63
C ITOL=11 is often useful for checking and comparing different
64
C routines. For this case, the user must supply the "exact"
65
C solution or a very accurate approximation (one with an error
66
C much less than TOL) through a common block,
67
C COMMON /SSLBLK/ SOLN( )
68
C If ITOL=11, iteration stops when the 2-norm of the difference
69
C between the iterative approximation and the user-supplied
70
C solution divided by the 2-norm of the user-supplied solution
71
C is less than TOL. Note that this requires the user to set up
72
C the "COMMON /SSLBLK/ SOLN(LENGTH)" in the calling routine.
73
C The routine with this declaration should be loaded before the
74
C stop test so that the correct length is used by the loader.
75
C This procedure is not standard Fortran and may not work
76
C correctly on your system (although it has worked on every
77
C system the authors have tried). If ITOL is not 11 then this
78
C common block is indeed standard Fortran.
80
C Convergence criterion, as described above. (Reset if IERR=4.)
82
C Maximum number of iterations.
84
C Number of iterations required to reach convergence, or
85
C ITMAX+1 if convergence criterion could not be achieved in
88
C Error estimate of error in final approximate solution, as
92
C IERR = 0 => All went well.
93
C IERR = 1 => Insufficient space allocated for WORK or IWORK.
94
C IERR = 2 => Method failed to converge in ITMAX steps.
95
C IERR = 3 => Error in user input.
96
C Check input values of N, ITOL.
97
C IERR = 4 => User error tolerance set too tight.
98
C Reset to 500*R1MACH(3). Iteration proceeded.
99
C IERR = 5 => Breakdown of the method detected.
100
C (r0,r) approximately 0.
101
C IERR = 6 => Stagnation of the method detected.
102
C (r0,v) approximately 0.
104
C Unit number on which to write the error at each iteration,
105
C if this is desired for monitoring convergence. If unit
106
C number is 0, no writing will occur.
107
C RWORK :WORK Real RWORK(LENW).
108
C Real array used for workspace.
110
C Length of the real workspace, RWORK. LENW >= 8*N.
111
C IWORK :WORK Integer IWORK(LENIW).
112
C Used to hold pointers into the RWORK array.
113
C Upon return the following locations of IWORK hold information
114
C which may be of use to the user:
115
C IWORK(9) Amount of Integer workspace actually used.
116
C IWORK(10) Amount of Real workspace actually used.
118
C Length of the integer workspace, IWORK. LENIW >= 10.
121
C This routine performs preconditioned BiConjugate gradient
122
C method on the Non-Symmetric positive definite linear system
123
C Ax=b. The preconditioner is M = DIAG(A), the diagonal of the
124
C matrix A. This is the simplest of preconditioners and
125
C vectorizes very well.
127
C The Sparse Linear Algebra Package (SLAP) utilizes two matrix
128
C data structures: 1) the SLAP Triad format or 2) the SLAP
129
C Column format. The user can hand this routine either of the
130
C of these data structures and SLAP will figure out which on
131
C is being used and act accordingly.
133
C =================== S L A P Triad format ===================
135
C This routine requires that the matrix A be stored in the
136
C SLAP Triad format. In this format only the non-zeros are
137
C stored. They may appear in *ANY* order. The user supplies
138
C three arrays of length NELT, where NELT is the number of
139
C non-zeros in the matrix: (IA(NELT), JA(NELT), A(NELT)). For
140
C each non-zero the user puts the row and column index of that
141
C matrix element in the IA and JA arrays. The value of the
142
C non-zero matrix element is placed in the corresponding
143
C location of the A array. This is an extremely easy data
144
C structure to generate. On the other hand it is not too
145
C efficient on vector computers for the iterative solution of
146
C linear systems. Hence, SLAP changes this input data
147
C structure to the SLAP Column format for the iteration (but
148
C does not change it back).
150
C Here is an example of the SLAP Triad storage format for a
151
C 5x5 Matrix. Recall that the entries may appear in any order.
153
C 5x5 Matrix SLAP Triad format for 5x5 matrix on left.
154
C 1 2 3 4 5 6 7 8 9 10 11
155
C |11 12 0 0 15| A: 51 12 11 33 15 53 55 22 35 44 21
156
C |21 22 0 0 0| IA: 5 1 1 3 1 5 5 2 3 4 2
157
C | 0 0 33 0 35| JA: 1 2 1 3 5 3 5 2 5 4 1
161
C =================== S L A P Column format ==================
163
C This routine requires that the matrix A be stored in the
164
C SLAP Column format. In this format the non-zeros are stored
165
C counting down columns (except for the diagonal entry, which
166
C must appear first in each "column") and are stored in the
167
C real array A. In other words, for each column in the matrix
168
C put the diagonal entry in A. Then put in the other non-zero
169
C elements going down the column (except the diagonal) in
170
C order. The IA array holds the row index for each non-zero.
171
C The JA array holds the offsets into the IA, A arrays for the
172
C beginning of each column. That is, IA(JA(ICOL)),
173
C A(JA(ICOL)) points to the beginning of the ICOL-th column in
174
C IA and A. IA(JA(ICOL+1)-1), A(JA(ICOL+1)-1) points to the
175
C end of the ICOL-th column. Note that we always have
176
C JA(N+1) = NELT+1, where N is the number of columns in the
177
C matrix and NELT is the number of non-zeros in the matrix.
179
C Here is an example of the SLAP Column storage format for a
180
C 5x5 Matrix (in the A and IA arrays '|' denotes the end of a
183
C 5x5 Matrix SLAP Column format for 5x5 matrix on left.
184
C 1 2 3 4 5 6 7 8 9 10 11
185
C |11 12 0 0 15| A: 11 21 51 | 22 12 | 33 53 | 44 | 55 15 35
186
C |21 22 0 0 0| IA: 1 2 5 | 2 1 | 3 5 | 4 | 5 1 3
187
C | 0 0 33 0 35| JA: 1 4 6 8 9 12
192
C The SLAP Triad format (IA, JA, A) is modified internally to
193
C be the SLAP Column format. See above.
196
C This routine will attempt to write to the Fortran logical output
197
C unit IUNIT, if IUNIT .ne. 0. Thus, the user must make sure that
198
C this logical unit is attached to a file or terminal before calling
199
C this routine with a non-zero value for IUNIT. This routine does
200
C not check for the validity of a non-zero IUNIT unit number.
202
C***SEE ALSO SCGS, SLUBCG
203
C***REFERENCES 1. P. Sonneveld, CGS, a fast Lanczos-type solver
204
C for nonsymmetric linear systems, Delft University
205
C of Technology Report 84-16, Department of Mathe-
206
C matics and Informatics, Delft, The Netherlands.
207
C 2. E. F. Kaasschieter, The solution of non-symmetric
208
C linear systems by biconjugate gradients or conjugate
209
C gradients squared, Delft University of Technology
210
C Report 86-21, Department of Mathematics and Informa-
211
C tics, Delft, The Netherlands.
212
C***ROUTINES CALLED SCGS, SCHKW, SS2Y, SSDI, SSDS, SSMV
213
C***REVISION HISTORY (YYMMDD)
214
C 871119 DATE WRITTEN
215
C 881213 Previous REVISION DATE
216
C 890915 Made changes requested at July 1989 CML Meeting. (MKS)
217
C 890921 Removed TeX from comments. (FNF)
218
C 890922 Numerous changes to prologue to make closer to SLATEC
220
C 890929 Numerous changes to reduce SP/DP differences. (FNF)
221
C 910411 Prologue converted to Version 4.0 format. (BAB)
222
C 920407 COMMON BLOCK renamed SSLBLK. (WRB)
223
C 920511 Added complete declaration section. (WRB)
224
C 920929 Corrected format of references. (FNF)
225
C 921113 Corrected C***CATEGORY line. (FNF)
226
C***END PROLOGUE SSDCGS
229
PARAMETER (LOCRB=1, LOCIB=11)
230
C .. Scalar Arguments ..
232
INTEGER IERR, ISYM, ITER, ITMAX, ITOL, IUNIT, LENIW, LENW, N, NELT
233
C .. Array Arguments ..
234
REAL A(N), B(N), RWORK(LENW), X(N)
235
INTEGER IA(NELT), IWORK(LENIW), JA(NELT)
236
C .. Local Scalars ..
237
INTEGER LOCDIN, LOCIW, LOCP, LOCQ, LOCR, LOCR0, LOCU, LOCV1,
239
C .. External Subroutines ..
240
EXTERNAL SCGS, SCHKW, SS2Y, SSDI, SSDS, SSMV
241
C***FIRST EXECUTABLE STATEMENT SSDCGS
244
IF( N.LT.1 .OR. NELT.LT.1 ) THEN
249
C Change the SLAP input matrix IA, JA, A to SLAP-Column format.
250
CALL SS2Y( N, NELT, IA, JA, A, ISYM )
252
C Set up the workspace.
265
C Check the workspace allocations.
266
CALL SCHKW( 'SSDCGS', LOCIW, LENIW, LOCW, LENW, IERR, ITER, ERR )
267
IF( IERR.NE.0 ) RETURN
273
C Compute the inverse of the diagonal of the matrix.
274
CALL SSDS(N, NELT, IA, JA, A, ISYM, RWORK(LOCDIN))
276
C Perform the Diagonally Scaled
277
C BiConjugate Gradient Squared algorithm.
278
CALL SCGS(N, B, X, NELT, IA, JA, A, ISYM, SSMV,
279
$ SSDI, ITOL, TOL, ITMAX, ITER, ERR, IERR, IUNIT,
280
$ RWORK(LOCR), RWORK(LOCR0), RWORK(LOCP),
281
$ RWORK(LOCQ), RWORK(LOCU), RWORK(LOCV1),
282
$ RWORK(LOCV2), RWORK(1), IWORK(1))
284
C------------- LAST LINE OF SSDCGS FOLLOWS ----------------------------