2
SUBROUTINE DSDCGN (N, B, X, NELT, IA, JA, A, ISYM, ITOL, TOL,
3
+ ITMAX, ITER, ERR, IERR, IUNIT, RWORK, LENW, IWORK, LENIW)
4
C***BEGIN PROLOGUE DSDCGN
5
C***PURPOSE Diagonally Scaled CG Sparse Ax=b Solver for Normal Eqn's.
6
C Routine to solve a general linear system Ax = b using
7
C diagonal scaling with the Conjugate Gradient method
8
C applied to the the normal equations, viz., AA'y = b,
10
C***LIBRARY SLATEC (SLAP)
11
C***CATEGORY D2A4, D2B4
12
C***TYPE DOUBLE PRECISION (SSDCGN-S, DSDCGN-D)
13
C***KEYWORDS ITERATIVE PRECONDITION, NON-SYMMETRIC LINEAR SYSTEM SOLVE,
15
C***AUTHOR Greenbaum, Anne, (Courant Institute)
16
C Seager, Mark K., (LLNL)
17
C Lawrence Livermore National Laboratory
19
C Livermore, CA 94550 (510) 423-3141
24
C INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX
25
C INTEGER ITER, IERR, IUNIT, LENW, IWORK(10), LENIW
26
C DOUBLE PRECISION B(N), X(N), A(NELT), TOL, ERR, RWORK(8*N)
28
C CALL DSDCGN(N, B, X, NELT, IA, JA, A, ISYM, ITOL, TOL,
29
C $ ITMAX, ITER, ERR, IERR, IUNIT, RWORK, LENW, IWORK, LENIW)
33
C Order of the Matrix.
34
C B :IN Double Precision B(N).
35
C Right-hand side vector.
36
C X :INOUT Double Precision X(N).
37
C On input X is your initial guess for solution vector.
38
C On output X is the final approximate solution.
40
C Number of Non-Zeros stored in A.
41
C IA :INOUT Integer IA(NELT).
42
C JA :INOUT Integer JA(NELT).
43
C A :INOUT Double Precision A(NELT).
44
C These arrays should hold the matrix A in either the SLAP
45
C Triad format or the SLAP Column format. See "Description",
46
C below. If the SLAP Triad format is chosen it is changed
47
C internally to the SLAP Column format.
49
C Flag to indicate symmetric storage format.
50
C If ISYM=0, all non-zero entries of the matrix are stored.
51
C If ISYM=1, the matrix is symmetric, and only the upper
52
C or lower triangle of the matrix is stored.
54
C Flag to indicate type of convergence criterion.
55
C If ITOL=1, iteration stops when the 2-norm of the residual
56
C divided by the 2-norm of the right-hand side is less than TOL.
57
C If ITOL=2, iteration stops when the 2-norm of M-inv times the
58
C residual divided by the 2-norm of M-inv times the right hand
59
C side is less than TOL, where M-inv is the inverse of the
61
C ITOL=11 is often useful for checking and comparing different
62
C routines. For this case, the user must supply the "exact"
63
C solution or a very accurate approximation (one with an error
64
C much less than TOL) through a common block,
65
C COMMON /DSLBLK/ SOLN( )
66
C If ITOL=11, iteration stops when the 2-norm of the difference
67
C between the iterative approximation and the user-supplied
68
C solution divided by the 2-norm of the user-supplied solution
69
C is less than TOL. Note that this requires the user to set up
70
C the "COMMON /DSLBLK/ SOLN(LENGTH)" in the calling routine.
71
C The routine with this declaration should be loaded before the
72
C stop test so that the correct length is used by the loader.
73
C This procedure is not standard Fortran and may not work
74
C correctly on your system (although it has worked on every
75
C system the authors have tried). If ITOL is not 11 then this
76
C common block is indeed standard Fortran.
77
C TOL :INOUT Double Precision.
78
C Convergence criterion, as described above. (Reset if IERR=4.)
80
C Maximum number of iterations.
82
C Number of iterations required to reach convergence, or
83
C ITMAX+1 if convergence criterion could not be achieved in
85
C ERR :OUT Double Precision.
86
C Error estimate of error in final approximate solution, as
90
C IERR = 0 => All went well.
91
C IERR = 1 => Insufficient space allocated for WORK or IWORK.
92
C IERR = 2 => Method failed to converge in ITMAX steps.
93
C IERR = 3 => Error in user input.
94
C Check input values of N, ITOL.
95
C IERR = 4 => User error tolerance set too tight.
96
C Reset to 500*D1MACH(3). Iteration proceeded.
97
C IERR = 5 => Preconditioning matrix, M, is not positive
98
C definite. (r,z) < 0.
99
C IERR = 6 => Matrix A is not positive definite. (p,Ap) < 0.
101
C Unit number on which to write the error at each iteration,
102
C if this is desired for monitoring convergence. If unit
103
C number is 0, no writing will occur.
104
C RWORK :WORK Double Precision RWORK(LENW).
105
C Double Precision array used for workspace.
107
C Length of the double precision workspace, RWORK.
109
C IWORK :WORK Integer IWORK(LENIW).
110
C Used to hold pointers into the RWORK array.
111
C Upon return the following locations of IWORK hold information
112
C which may be of use to the user:
113
C IWORK(9) Amount of Integer workspace actually used.
114
C IWORK(10) Amount of Double Precision workspace actually used.
116
C Length of the integer workspace, IWORK. LENIW >= 10.
119
C This routine is simply a driver for the DCGN routine. It
120
C calls the DSD2S routine to set up the preconditioning and
121
C then calls DCGN with the appropriate MATVEC and MSOLVE
124
C The Sparse Linear Algebra Package (SLAP) utilizes two matrix
125
C data structures: 1) the SLAP Triad format or 2) the SLAP
126
C Column format. The user can hand this routine either of the
127
C of these data structures and SLAP will figure out which on
128
C is being used and act accordingly.
130
C =================== S L A P Triad format ===================
132
C This routine requires that the matrix A be stored in the
133
C SLAP Triad format. In this format only the non-zeros are
134
C stored. They may appear in *ANY* order. The user supplies
135
C three arrays of length NELT, where NELT is the number of
136
C non-zeros in the matrix: (IA(NELT), JA(NELT), A(NELT)). For
137
C each non-zero the user puts the row and column index of that
138
C matrix element in the IA and JA arrays. The value of the
139
C non-zero matrix element is placed in the corresponding
140
C location of the A array. This is an extremely easy data
141
C structure to generate. On the other hand it is not too
142
C efficient on vector computers for the iterative solution of
143
C linear systems. Hence, SLAP changes this input data
144
C structure to the SLAP Column format for the iteration (but
145
C does not change it back).
147
C Here is an example of the SLAP Triad storage format for a
148
C 5x5 Matrix. Recall that the entries may appear in any order.
150
C 5x5 Matrix SLAP Triad format for 5x5 matrix on left.
151
C 1 2 3 4 5 6 7 8 9 10 11
152
C |11 12 0 0 15| A: 51 12 11 33 15 53 55 22 35 44 21
153
C |21 22 0 0 0| IA: 5 1 1 3 1 5 5 2 3 4 2
154
C | 0 0 33 0 35| JA: 1 2 1 3 5 3 5 2 5 4 1
158
C =================== S L A P Column format ==================
160
C This routine requires that the matrix A be stored in the
161
C SLAP Column format. In this format the non-zeros are stored
162
C counting down columns (except for the diagonal entry, which
163
C must appear first in each "column") and are stored in the
164
C double precision array A. In other words, for each column
165
C in the matrix put the diagonal entry in A. Then put in the
166
C other non-zero elements going down the column (except the
167
C diagonal) in order. The IA array holds the row index for
168
C each non-zero. The JA array holds the offsets into the IA,
169
C A arrays for the beginning of each column. That is,
170
C IA(JA(ICOL)), A(JA(ICOL)) points to the beginning of the
171
C ICOL-th column in IA and A. IA(JA(ICOL+1)-1),
172
C A(JA(ICOL+1)-1) points to the end of the ICOL-th column.
173
C Note that we always have JA(N+1) = NELT+1, where N is the
174
C number of columns in the matrix and NELT is the number of
175
C non-zeros in the matrix.
177
C Here is an example of the SLAP Column storage format for a
178
C 5x5 Matrix (in the A and IA arrays '|' denotes the end of a
181
C 5x5 Matrix SLAP Column format for 5x5 matrix on left.
182
C 1 2 3 4 5 6 7 8 9 10 11
183
C |11 12 0 0 15| A: 11 21 51 | 22 12 | 33 53 | 44 | 55 15 35
184
C |21 22 0 0 0| IA: 1 2 5 | 2 1 | 3 5 | 4 | 5 1 3
185
C | 0 0 33 0 35| JA: 1 4 6 8 9 12
190
C The SLAP Triad format (IA, JA, A) is modified internally to be
191
C the SLAP Column format. See above.
194
C This routine will attempt to write to the Fortran logical output
195
C unit IUNIT, if IUNIT .ne. 0. Thus, the user must make sure that
196
C this logical unit is attached to a file or terminal before calling
197
C this routine with a non-zero value for IUNIT. This routine does
198
C not check for the validity of a non-zero IUNIT unit number.
200
C***SEE ALSO DCGN, DSD2S, DSMV, DSMTV, DSDI
201
C***REFERENCES (NONE)
202
C***ROUTINES CALLED DCGN, DCHKW, DS2Y, DSD2S, DSDI, DSMTV, DSMV
203
C***REVISION HISTORY (YYMMDD)
204
C 890404 DATE WRITTEN
205
C 890404 Previous REVISION DATE
206
C 890915 Made changes requested at July 1989 CML Meeting. (MKS)
207
C 890921 Removed TeX from comments. (FNF)
208
C 890922 Numerous changes to prologue to make closer to SLATEC
210
C 890929 Numerous changes to reduce SP/DP differences. (FNF)
211
C 910411 Prologue converted to Version 4.0 format. (BAB)
212
C 920407 COMMON BLOCK renamed DSLBLK. (WRB)
213
C 920511 Added complete declaration section. (WRB)
214
C 921113 Corrected C***CATEGORY line. (FNF)
215
C***END PROLOGUE DSDCGN
218
PARAMETER (LOCRB=1, LOCIB=11)
219
C .. Scalar Arguments ..
220
DOUBLE PRECISION ERR, TOL
221
INTEGER IERR, ISYM, ITER, ITMAX, ITOL, IUNIT, LENIW, LENW, N, NELT
222
C .. Array Arguments ..
223
DOUBLE PRECISION A(NELT), B(N), RWORK(LENW), X(N)
224
INTEGER IA(NELT), IWORK(LENIW), JA(NELT)
225
C .. Local Scalars ..
226
INTEGER LOCATD, LOCATP, LOCATZ, LOCD, LOCDZ, LOCIW, LOCP, LOCR,
228
C .. External Subroutines ..
229
EXTERNAL DCGN, DCHKW, DS2Y, DSD2S, DSDI, DSMTV, DSMV
230
C***FIRST EXECUTABLE STATEMENT DSDCGN
233
IF( N.LT.1 .OR. NELT.LT.1 ) THEN
238
C Modify the SLAP matrix data structure to YSMP-Column.
239
CALL DS2Y( N, NELT, IA, JA, A, ISYM )
241
C Set up the work arrays.
254
C Check the workspace allocations.
255
CALL DCHKW( 'DSDCGN', LOCIW, LENIW, LOCW, LENW, IERR, ITER, ERR )
256
IF( IERR.NE.0 ) RETURN
262
C Compute the inverse of the diagonal of AA'. This will be
263
C used as the preconditioner.
264
CALL DSD2S(N, NELT, IA, JA, A, ISYM, RWORK(1))
266
C Perform Conjugate Gradient algorithm on the normal equations.
267
CALL DCGN( N, B, X, NELT, IA, JA, A, ISYM, DSMV, DSMTV, DSDI,
268
$ ITOL, TOL, ITMAX, ITER, ERR, IERR, IUNIT, RWORK(LOCR),
269
$ RWORK(LOCZ), RWORK(LOCP), RWORK(LOCATP), RWORK(LOCATZ),
270
$ RWORK(LOCDZ), RWORK(LOCATD), RWORK, IWORK )
272
IF( ITER.GT.ITMAX ) IERR = 2
274
C------------- LAST LINE OF DSDCGN FOLLOWS ----------------------------
added
removed
deps/openlibm/slatec/dsdcgn.f
