2
SUBROUTINE <_c>GMRESREVCOM(N, B, X, RESTRT, WORK, LDW, WORK2,
3
$ LDW2, ITER, RESID, INFO, NDX1, NDX2, SCLR1,
6
* -- Iterative template routine --
7
* Univ. of Tennessee and Oak Ridge National Laboratory
9
* Details of this algorithm are described in "Templates for the
10
* Solution of Linear Systems: Building Blocks for Iterative
11
* Methods", Barrett, Berry, Chan, Demmel, Donato, Dongarra,
12
* EiITERkhout, Pozo, Romine, and van der Vorst, SIAM Publications,
13
* 1993. (ftp netlib2.cs.utk.edu; cd linalg; get templates.ps).
15
* .. Scalar Arguments ..
16
INTEGER N, RESTRT, LDW, LDW2, ITER, INFO
17
<rt=real,double precision,real,double precision> RESID
22
* .. Array Arguments ..
23
<_t> B( * ), X( * ), WORK( LDW,* ), WORK2( LDW2,* )
29
* GMRES solves the linear system Ax = b using the
30
* Generalized Minimal Residual iterative method with preconditioning.
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; on exit, the iterated solution.
46
* RESTRT (input) INTEGER
47
* Restart parameter, .ls. = N. This parameter controls the amount
48
* of memory required for matrix WORK2.
50
* WORK (workspace) DOUBLE PRECISION array, dimension (LDW,6+restrt).
51
* Note that if the initial guess is the zero vector, then
52
* storing the initial residual is not necessary.
55
* The leading dimension of the array WORK. LDW .gt. = max(1,N).
57
* WORK2 (workspace) DOUBLE PRECISION array, dimension (LDW2,2*RESTRT+2).
58
* This workspace is used for constructing and storing the
59
* upper Hessenberg matrix. The two extra columns are used to
60
* store the Givens rotation matrices.
62
* LDW2 (input) INTEGER
63
* The leading dimension of the array WORK2.
64
* LDW2 .gt. = max(2,RESTRT+1).
66
* ITER (input/output) INTEGER
67
* On input, the maximum iterations to be performed.
68
* On output, actual number of iterations performed.
70
* RESID (input/output) DOUBLE PRECISION
71
* On input, the allowable error tolerance.
72
* On ouput, the norm of the residual vector if solution
73
* approximated to tolerance, otherwise reset to input
76
* INFO (output) INTEGER
77
* = 0: successful exit
78
* = 1: maximum number of iterations performed;
79
* convergence not achieved.
80
* -5: Erroneous NDX1/NDX2 in INIT call.
83
* NDX1 (input/output) INTEGER.
84
* NDX2 On entry in INIT call contain indices required by interface
85
* level for stopping test.
86
* All other times, used as output, to indicate indices into
87
* WORK[] for the MATVEC, PSOLVE done by the interface level.
89
* SCLR1 (output) DOUBLE PRECISION.
90
* SCLR2 Used to pass the scalars used in MATVEC. Scalars are reqd because
91
* original routines use dgemv.
93
* IJOB (input/output) INTEGER.
94
* Used to communicate job code between the two levels.
96
* ============================================================
100
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
102
PARAMETER ( OFSET = 1000 )
104
* .. Local Scalars ..
105
INTEGER I, MAXIT, AV, GIV, H, R, S, V, W, Y,
107
<_t> <xdot=sdot,ddot,cdotc,zdotc>
109
<rt> BNRM2, RNORM, TOL,
110
$ <rc=s,d,sc,dz>NRM2,
114
* indicates where to resume from. Only valid when IJOB = 2!
121
* .. External Routines ..
122
EXTERNAL <_c>AXPY, <_c>COPY, <xdot>, <rc>NRM2, <_c>SCAL
124
* .. Executable Statements ..
126
* Entry point, so test IJOB
127
IF (IJOB .eq. 1) THEN
129
ELSEIF (IJOB .eq. 2) THEN
130
* here we do resumption handling
131
IF (RLBL .eq. 2) GOTO 2
132
IF (RLBL .eq. 3) GOTO 3
133
IF (RLBL .eq. 4) GOTO 4
134
IF (RLBL .eq. 5) GOTO 5
135
IF (RLBL .eq. 6) GOTO 6
136
IF (RLBL .eq. 7) GOTO 7
137
* if neither of these, then error
151
* Alias workspace columns.
163
* Check if caller will need indexing info.
165
IF( NDX1.NE.-1 ) THEN
167
NEED1 = ((R - 1) * LDW) + 1
168
ELSEIF( NDX1.EQ.2 ) THEN
169
NEED1 = ((S - 1) * LDW) + 1
170
ELSEIF( NDX1.EQ.3 ) THEN
171
NEED1 = ((W - 1) * LDW) + 1
172
ELSEIF( NDX1.EQ.4 ) THEN
173
NEED1 = ((Y - 1) * LDW) + 1
174
ELSEIF( NDX1.EQ.5 ) THEN
175
NEED1 = ((AV - 1) * LDW) + 1
176
ELSEIF( NDX1.EQ.6 ) THEN
177
NEED1 = ((V - 1) * LDW) + 1
178
ELSEIF( ( NDX1.GT.V*OFSET ) .AND.
179
$ ( NDX1.LE.V*OFSET+RESTRT ) ) THEN
180
NEED1 = ((NDX1-V*OFSET - 1) * LDW) + 1
181
ELSEIF( ( NDX1.GT.GIV*OFSET ) .AND.
182
$ ( NDX1.LE.GIV*OFSET+RESTRT ) ) THEN
183
NEED1 = ((NDX1-GIV*OFSET - 1) * LDW) + 1
193
IF( NDX2.NE.-1 ) THEN
195
NEED2 = ((R - 1) * LDW) + 1
196
ELSEIF( NDX2.EQ.2 ) THEN
197
NEED2 = ((S - 1) * LDW) + 1
198
ELSEIF( NDX2.EQ.3 ) THEN
199
NEED2 = ((W - 1) * LDW) + 1
200
ELSEIF( NDX2.EQ.4 ) THEN
201
NEED2 = ((Y - 1) * LDW) + 1
202
ELSEIF( NDX2.EQ.5 ) THEN
203
NEED2 = ((AV - 1) * LDW) + 1
204
ELSEIF( NDX2.EQ.6 ) THEN
205
NEED2 = ((V - 1) * LDW) + 1
206
ELSEIF( ( NDX2.GT.V*OFSET ) .AND.
207
$ ( NDX2.LE.V*OFSET+RESTRT ) ) THEN
208
NEED2 = ((NDX2-V*OFSET - 1) * LDW) + 1
209
ELSEIF( ( NDX2.GT.GIV*OFSET ) .AND.
210
$ ( NDX2.LE.GIV*OFSET+RESTRT ) ) THEN
211
NEED2 = ((NDX2-GIV*OFSET - 1) * LDW) + 1
221
* Set initial residual.
223
CALL <_c>COPY( N, B, 1, WORK(1,R), 1 )
224
IF ( <rc>NRM2( N, X, 1 ).NE.ZERO ) THEN
225
*********CALL MATVEC( -ONE, X, ONE, WORK(1,R) )
226
* Note: using X directly
230
NDX2 = ((R - 1) * LDW) + 1
232
* Prepare for resumption & return
241
IF ( <rc>NRM2( N, WORK(1,R), 1 ).LT.TOL ) GO TO 200
243
BNRM2 = <rc>NRM2( N, B, 1 )
244
IF ( BNRM2.EQ.ZERO ) BNRM2 = ONE
251
* Construct the first column of V, and initialize S to the
252
* elementary vector E1 scaled by RNORM.
254
*********CALL PSOLVE( WORK( 1,V ), WORK( 1,R ) )
256
NDX1 = ((V - 1) * LDW) + 1
257
NDX2 = ((R - 1) * LDW) + 1
258
* Prepare for return & return
267
RNORM = <rc>NRM2( N, WORK( 1,V ), 1 )
269
CALL <_c>SCAL( N, toz, WORK( 1,V ), 1 )
270
CALL <_c>ELEMVEC( 1, N, RNORM, WORK( 1,S ) )
272
* DO 50 I = 1, RESTRT
274
49 if (i.gt.restrt) go to 50
275
************CALL MATVEC( ONE, WORK( 1,V+I-1 ), ZERO, WORK( 1,AV ) )
277
NDX1 = ((V+I-1 - 1) * LDW) + 1
278
NDX2 = ((AV - 1) * LDW) + 1
279
* Prepare for return & return
290
*********CALL PSOLVE( WORK( 1,W ), WORK( 1,AV ) )
292
NDX1 = ((W - 1) * LDW) + 1
293
NDX2 = ((AV - 1) * LDW) + 1
294
* Prepare for return & return
303
* Construct I-th column of H so that it is orthnormal to
304
* the previous I-1 columns.
306
CALL <_c>ORTHOH( I, N, WORK2( 1,I+H-1 ), WORK( 1,V ), LDW,
311
* Apply Givens rotations to the I-th column of H. This
312
* effectively reduces the Hessenberg matrix to upper
313
* triangular form during the RESTRT iterations.
315
$ CALL <_c>APPLYGIVENS(I, WORK2( 1,I+H-1 ), WORK2( 1,GIV ),
318
* Approxiyate residual norm. Check tolerance. If okay, compute
319
* final approximation vector X and quit.
321
RESID = <rc>APPROXRES( I, WORK2( 1,I+H-1 ), WORK( 1,S ),
322
$ WORK2( 1,GIV ), LDW2 ) / BNRM2
323
IF ( RESID.LE.TOL ) THEN
324
CALL <_c>UPDATE(I, N, X, WORK2( 1,H ), LDW2,
325
$ WORK(1,Y), WORK(1,S), WORK( 1,V ), LDW)
333
* Compute current solution vector X.
335
CALL <_c>UPDATE(RESTRT, N, X, WORK2( 1,H ), LDW2,
336
$ WORK(1,Y), WORK( 1,S ), WORK( 1,V ), LDW )
338
* Compute residual vector R, find norm,
339
* then check for tolerance.
341
CALL <_c>COPY( N, B, 1, WORK( 1,R ), 1 )
342
*********CALL MATVEC( -ONE, X, ONE, WORK( 1,R ) )
345
NDX2 = ((R - 1) * LDW) + 1
346
* Prepare for return & return
357
WORK( I+1,S ) = <rc>NRM2( N, WORK( 1,R ), 1 )
359
*********RESID = WORK( I+1,S ) / BNRM2
360
*********IF ( RESID.LE.TOL ) GO TO 200
364
* Prepare for resumption & return
372
IF( INFO.EQ.1 ) GO TO 200
374
IF ( ITER.EQ.MAXIT ) THEN
391
* Iteration successful; return.
401
END SUBROUTINE <_c>GMRESREVCOM
403
* =========================================================
404
SUBROUTINE <_c>ORTHOH( I, N, H, V, LDV, W )
407
<_t> H( * ), W( * ), V( LDV,* )
409
* Construct the I-th column of the upper Hessenberg matrix H
410
* using the Gram-Schmidt process on V and W.
413
<rt=real,double precision,real,double precision>
414
$ <rc=s,d,sc,dz>NRM2, ONE
415
PARAMETER ( ONE = 1.0D+0 )
416
<_t> <xdot=sdot,ddot,cdotc,zdotc>
417
EXTERNAL <_c>AXPY, <_c>COPY, <xdot>, <rc>NRM2, <_c>SCAL
420
H( K ) = <xdot>( N, V( 1,K ), 1, W, 1 )
421
CALL <_c>AXPY( N, -H( K ), V( 1,K ), 1, W, 1 )
423
H( I+1 ) = <rc>NRM2( N, W, 1 )
424
CALL <_c>COPY( N, W, 1, V( 1,I+1 ), 1 )
425
CALL <_c>SCAL( N, ONE / H( I+1 ), V( 1,I+1 ), 1 )
429
END SUBROUTINE <_c>ORTHOH
430
* =========================================================
431
SUBROUTINE <_c>APPLYGIVENS( I, H, GIVENS, LDG )
434
<_t> H( * ), GIVENS( LDG,* )
436
* This routine applies a sequence of I-1 Givens rotations to
437
* the I-th column of H. The Givens parameters are stored, so that
438
* the first I-2 Givens rotation matrices are known. The I-1st
439
* Givens rotation is computed using BLAS 1 routine DROTG. Each
440
* rotation is applied to the 2x1 vector [H( J ), H( J+1 )]',
441
* which results in H( J+1 ) = 0.
444
* DOUBLE PRECISION TEMP
447
* .. Executable Statements ..
449
* Construct I-1st rotation matrix.
451
* CALL <_c>ROTG( H( I ), H( I+1 ), GIVENS( I,1 ), GIVENS( I,2 ) )
452
* CALL <_c>GETGIV( H( I ), H( I+1 ), GIVENS( I,1 ), GIVENS( I,2 ) )
454
* Apply 1,...,I-1st rotation matrices to the I-th column of H.
457
CALL <_c>ROTVEC(H( J ), H( J+1 ), GIVENS( J,1 ), GIVENS( J,2 ))
458
* TEMP = GIVENS( J,1 ) * H( J ) + GIVENS( J,2 ) * H( J+1 )
459
* H( J+1 ) = -GIVENS( J,2 ) * H( J ) + GIVENS( J,1 ) * H( J+1 )
462
call <_c>getgiv( H( I ), H( I+1 ), GIVENS( I,1 ), GIVENS( I,2 ) )
463
call <_c>rotvec( H( I ), H( I+1 ), GIVENS( I,1 ), GIVENS( I,2 ) )
467
END SUBROUTINE <_c>APPLYGIVENS
469
* ===============================================================
470
<rt=real,double precision,real,double precision>
471
$ FUNCTION <rc=s,d,sc,dz>APPROXRES( I, H, S, GIVENS, LDG )
474
<_t> H( * ), S( * ), GIVENS( LDG,* )
476
* This func allows the user to approximate the residual
477
* using an updating scheme involving Givens rotations. The
478
* rotation matrix is formed using [H( I ),H( I+1 )]' with the
479
* intent of zeroing H( I+1 ), but here is applied to the 2x1
480
* vector [S(I), S(I+1)]'.
485
* .. Executable Statements ..
487
* CALL <_c>ROTG( H( I ), H( I+1 ), GIVENS( I,1 ), GIVENS( I,2 ) )
488
* CALL <_c>GETGIV( H( I ), H( I+1 ), GIVENS( I,1 ), GIVENS( I,2 ) )
489
CALL <_c>ROTVEC( S( I ), S( I+1 ), GIVENS( I,1 ), GIVENS( I,2 ) )
491
<rc>APPROXRES = ABS( S( I+1 ) )
495
END FUNCTION <rc>APPROXRES
496
* ===============================================================
497
SUBROUTINE <_c>UPDATE( I, N, X, H, LDH, Y, S, V, LDV )
499
INTEGER N, I, J, LDH, LDV
500
<_t> X( * ), Y( * ), S( * ), H( LDH,* ), V( LDV,* )
501
EXTERNAL <_c>AXPY, <_c>COPY, <_c>TRSV
503
* Solve H*y = s for upper triangualar H.
505
CALL <_c>COPY( I, S, 1, Y, 1 )
506
CALL <_c>TRSV( 'UPPER', 'NOTRANS', 'NONUNIT', I, H, LDH, Y, 1 )
508
* Compute current solution vector X.
511
CALL <_c>AXPY( N, Y( J ), V( 1,J ), 1, X, 1 )
516
END SUBROUTINE <_c>UPDATE
518
* ===============================================================
519
SUBROUTINE <_c>GETGIV( A, B, C, S )
521
<_t> A, B, C, S, TEMP, ZERO, ONE
526
IF ( ABS( B ).EQ.ZERO ) THEN
529
ELSE IF ( ABS( B ).GT.ABS( A ) ) THEN
531
S = ONE / SQRT( ONE + abs(TEMP)**2 )
533
* S = b / SQRT( abs(a)**2 + abs(b)**2 )
534
* C = -a / SQRT( abs(a)**2 + abs(b)**2 )
537
C = ONE / SQRT( ONE + abs(TEMP)**2 )
539
* S = -b / SQRT( abs(a)**2 + abs(b)**2 )
540
* C = a / SQRT( abs(a)**2 + abs(b)**2 )
545
END SUBROUTINE <_c>GETGIV
547
* ================================================================
548
SUBROUTINE <_c>ROTVEC( X, Y, C, S )
550
<_t> X, Y, C, S, TEMP
553
TEMP = <co= , ,conjg,conjg>(C) * X - <co>(S) * Y
559
END SUBROUTINE <_c>ROTVEC
561
* ===============================================================
562
SUBROUTINE <_c>ELEMVEC( I, N, ALPHA, E )
564
* Construct the I-th elementary vector E, scaled by ALPHA.
570
<rt=real,double precision,real,double precision> ZERO
571
PARAMETER ( ZERO = 0.0D+0 )
580
END SUBROUTINE <_c>ELEMVEC
added
removed
Lib/linalg/iterative/GMRESREVCOM.f.src
