2
SUBROUTINE <_c>BICGSTABREVCOM(N, B, X, WORK, LDW, ITER, RESID,
3
$ INFO,NDX1, NDX2, SCLR1, SCLR2, IJOB)
5
* -- Iterative template routine --
6
* Univ. of Tennessee and Oak Ridge National Laboratory
8
* Details of this algorithm are described in "Templates for the
9
* Solution of Linear Systems: Building Blocks for Iterative
10
* Methods", Barrett, Berry, Chan, Demmel, Donato, Dongarra,
11
* Eijkhout, Pozo, Romine, and van der Vorst, SIAM Publications,
12
* 1993. (ftp netlib2.cs.utk.edu; cd linalg; get templates.ps).
14
* .. Scalar Arguments ..
15
INTEGER N, LDW, ITER, INFO
16
<rt=real,double precision,real,double precision> RESID
21
* .. Array Arguments ..
22
<_t> X( * ), B( * ), WORK( LDW,* )
28
* BICGSTAB solves the linear system A*x = b using the
29
* BiConjugate Gradient Stabilized iterative method with
36
* On entry, the dimension of the matrix.
39
* B (input) DOUBLE PRECISION array, dimension N.
40
* On entry, right hand side vector B.
43
* X (input/output) DOUBLE PRECISION array, dimension N.
44
* On input, the initial guess. This is commonly set to
46
* On exit, if INFO = 0, the iterated approximate solution.
48
* WORK (workspace) DOUBLE PRECISION array, dimension (LDW,7)
49
* Workspace for residual, direction vector, etc.
50
* Note that vectors R and S shared the same workspace.
53
* The leading dimension of the array WORK. LDW .gt. = max(1,N).
55
* ITER (input/output) INTEGER
56
* On input, the maximum iterations to be performed.
57
* On output, actual number of iterations performed.
59
* RESID (input/output) DOUBLE PRECISION
60
* On input, the allowable convergence measure for
61
* norm( b - A*x ) / norm( b ).
62
* On output, the final value of this measure.
64
* INFO (output) INTEGER
66
* = 0: Successful exit. Iterated approximate solution returned.
68
* .gt. 0: Convergence to tolerance not achieved. This will be
69
* set to the number of iterations performed.
71
* .ls. 0: Illegal input parameter, or breakdown occurred
76
* -1: matrix dimension N .ls. 0
78
* -3: Maximum number of iterations ITER .ls. = 0.
79
* -5: Erroneous NDX1/NDX2 in INIT call.
82
* BREAKDOWN: If parameters RHO or OMEGA become smaller
83
* than some tolerance, the program will terminate.
84
* Here we check against tolerance BREAKTOL.
86
* -10: RHO .ls. BREAKTOL: RHO and RTLD have become
88
* -11: OMEGA .ls. BREAKTOL: S and T have become
89
* orthogonal relative to T'*T.
91
* BREAKTOL is set in func GETBREAK.
93
* NDX1 (input/output) INTEGER.
94
* NDX2 On entry in INIT call contain indices required by interface
95
* level for stopping test.
96
* All other times, used as output, to indicate indices into
97
* WORK[] for the MATVEC, PSOLVE done by the interface level.
99
* SCLR1 (output) DOUBLE PRECISION.
100
* SCLR2 Used to pass the scalars used in MATVEC. Scalars are reqd because
101
* original routines use dgemv.
103
* IJOB (input/output) INTEGER.
104
* Used to communicate job code between the two levels.
106
* BLAS CALLS: DAXPY, DCOPY, DDOT, DNRM2, DSCAL
107
* ==============================================================
110
DOUBLE PRECISION ZERO, ONE
111
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
113
* .. Local Scalars ..
114
INTEGER R, RTLD, P, PHAT, V, S, SHAT, T, MAXIT,
117
$ RHOTOL, OMEGATOL, <sdsd=s,d,s,d>GETBREAK,
119
<_t> ALPHA, BETA, RHO, RHO1, OMEGA,
120
$ <xdot=sdot,ddot,cdotc,zdotc>
121
* indicates where to resume from. Only valid when IJOB = 2!
127
* .. External Funcs ..
128
EXTERNAL <sdsd>GETBREAK, <_c>AXPY, <_c>COPY,
129
$ <xdot>, <rc>NRM2, <_c>SCAL
131
* .. Intrinsic Funcs ..
134
* .. Executable Statements ..
136
* Entry point, so test IJOB
137
IF (IJOB .eq. 1) THEN
139
ELSEIF (IJOB .eq. 2) THEN
140
* here we do resumption handling
141
IF (RLBL .eq. 2) GOTO 2
142
IF (RLBL .eq. 3) GOTO 3
143
IF (RLBL .eq. 4) GOTO 4
144
IF (RLBL .eq. 5) GOTO 5
145
IF (RLBL .eq. 6) GOTO 6
146
IF (RLBL .eq. 7) GOTO 7
147
* if neither of these, then error
161
* Alias workspace columns.
172
* Check if caller will need indexing info.
174
IF( NDX1.NE.-1 ) THEN
176
NEED1 = ((R - 1) * LDW) + 1
177
ELSEIF( NDX1.EQ.2 ) THEN
178
NEED1 = ((RTLD - 1) * LDW) + 1
179
ELSEIF( NDX1.EQ.3 ) THEN
180
NEED1 = ((P - 1) * LDW) + 1
181
ELSEIF( NDX1.EQ.4 ) THEN
182
NEED1 = ((V - 1) * LDW) + 1
183
ELSEIF( NDX1.EQ.5 ) THEN
184
NEED1 = ((T - 1) * LDW) + 1
185
ELSEIF( NDX1.EQ.6 ) THEN
186
NEED1 = ((PHAT - 1) * LDW) + 1
187
ELSEIF( NDX1.EQ.7 ) THEN
188
NEED1 = ((SHAT - 1) * LDW) + 1
189
ELSEIF( NDX1.EQ.8 ) THEN
190
NEED1 = ((S - 1) * LDW) + 1
200
IF( NDX2.NE.-1 ) THEN
202
NEED2 = ((R - 1) * LDW) + 1
203
ELSEIF( NDX2.EQ.2 ) THEN
204
NEED2 = ((RTLD - 1) * LDW) + 1
205
ELSEIF( NDX2.EQ.3 ) THEN
206
NEED2 = ((P - 1) * LDW) + 1
207
ELSEIF( NDX2.EQ.4 ) THEN
208
NEED2 = ((V - 1) * LDW) + 1
209
ELSEIF( NDX2.EQ.5 ) THEN
210
NEED2 = ((T - 1) * LDW) + 1
211
ELSEIF( NDX2.EQ.6 ) THEN
212
NEED2 = ((PHAT - 1) * LDW) + 1
213
ELSEIF( NDX2.EQ.7 ) THEN
214
NEED2 = ((SHAT - 1) * LDW) + 1
215
ELSEIF( NDX2.EQ.8 ) THEN
216
NEED2 = ((S - 1) * LDW) + 1
226
* Set parameter tolerances.
228
RHOTOL = <sdsd>GETBREAK()
229
OMEGATOL = <sdsd>GETBREAK()
231
* Set initial residual.
233
CALL <_c>COPY( N, B, 1, WORK(1,R), 1 )
234
IF ( <rc>NRM2( N, X, 1 ).NE.ZERO ) THEN
235
*********CALL <_c>MATVEC( -ONE, X, ONE, WORK(1,R) )
236
* Note: using RTLD[] as temp. storage.
237
*********CALL <_c>COPY(N, X, 1, WORK(1,RTLD), 1)
241
NDX2 = ((R - 1) * LDW) + 1
243
* Prepare for resumption & return
252
IF ( <rc>NRM2( N, WORK(1,R), 1 ).LE.TOL ) GO TO 30
254
CALL <_c>COPY( N, WORK(1,R), 1, WORK(1,RTLD), 1 )
256
BNRM2 = <rc>NRM2( N, B, 1 )
257
IF ( BNRM2 .EQ. ZERO ) BNRM2 = ONE
263
* Perform BiConjugate Gradient Stabilized iteration.
267
RHO = <xdot>( N, WORK(1,RTLD), 1, WORK(1,R), 1 )
268
IF ( ABS( RHO ).LT.RHOTOL ) GO TO 25
272
IF ( ITER.GT.1 ) THEN
273
BETA = ( RHO / RHO1 ) * ( ALPHA / OMEGA )
274
CALL <_c>AXPY( N, -OMEGA, WORK(1,V), 1, WORK(1,P), 1 )
275
CALL <_c>SCAL( N, BETA, WORK(1,P), 1 )
276
CALL <_c>AXPY( N, ONE, WORK(1,R), 1, WORK(1,P), 1 )
278
CALL <_c>COPY( N, WORK(1,R), 1, WORK(1,P), 1 )
281
* Compute direction adjusting vector PHAT and scalar ALPHA.
283
*********CALL PSOLVE( WORK(1,PHAT), WORK(1,P) )
285
NDX1 = ((PHAT - 1) * LDW) + 1
286
NDX2 = ((P - 1) * LDW) + 1
287
* Prepare for return & return
296
*********CALL MATVEC( ONE, WORK(1,PHAT), ZERO, WORK(1,V) )
298
NDX1 = ((PHAT - 1) * LDW) + 1
299
NDX2 = ((V - 1) * LDW) + 1
300
* Prepare for return & return
311
ALPHA = RHO / <xdot>( N, WORK(1,RTLD), 1, WORK(1,V), 1 )
313
* Early check for tolerance.
315
CALL <_c>AXPY( N, -ALPHA, WORK(1,V), 1, WORK(1,R), 1 )
316
CALL <_c>COPY( N, WORK(1,R), 1, WORK(1,S), 1 )
317
IF ( <rc>NRM2( N, WORK(1,S), 1 ).LE.TOL ) THEN
318
CALL <_c>AXPY( N, ALPHA, WORK(1,PHAT), 1, X, 1 )
319
RESID = <rc>NRM2( N, WORK(1,S), 1 ) / BNRM2
323
* Compute stabilizer vector SHAT and scalar OMEGA.
325
************CALL PSOLVE( WORK(1,SHAT), WORK(1,S) )
327
NDX1 = ((SHAT - 1) * LDW) + 1
328
NDX2 = ((S - 1) * LDW) + 1
329
* Prepare for return & return
338
************CALL MATVEC( ONE, WORK(1,SHAT), ZERO, WORK(1,T) )
340
NDX1 = ((SHAT - 1) * LDW) + 1
341
NDX2 = ((T - 1) * LDW) + 1
342
* Prepare for return & return
353
OMEGA = <xdot>( N, WORK(1,T), 1, WORK(1,S), 1 ) /
354
$ <xdot>( N, WORK(1,T), 1, WORK(1,T), 1 )
356
* Compute new solution approximation vector X.
358
CALL <_c>AXPY( N, ALPHA, WORK(1,PHAT), 1, X, 1 )
359
CALL <_c>AXPY( N, OMEGA, WORK(1,SHAT), 1, X, 1 )
361
* Compute residual R, check for tolerance.
363
CALL <_c>AXPY( N, -OMEGA, WORK(1,T), 1, WORK(1,R), 1 )
365
************RESID = DNRM2( N, WORK(1,R), 1 ) / BNRM2
366
************IF ( RESID.LE.TOL ) GO TO 30
370
* Prepare for resumption & return
378
IF( INFO.EQ.1 ) GO TO 30
380
IF ( ITER.EQ.MAXIT ) THEN
387
IF ( ABS( OMEGA ).LT.OMEGATOL ) THEN
404
* Set breakdown flag.
406
IF ( ABS( RHO ).LT.RHOTOL ) THEN
408
ELSE IF ( ABS( OMEGA ).LT.OMEGATOL ) THEN
417
* Iteration successful; return.
424
* End of BICGSTABREVCOM
426
END SUBROUTINE <_c>BICGSTABREVCOM
added
removed
Lib/linalg/iterative/BiCGSTABREVCOM.f.src
