← Back to branch summary

~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

Viewing changes to Lib/linalg/iterative/BiCGREVCOM.f.src

  • Committer: Bazaar Package Importer
  • Author(s): Daniel T. Chen (new)
  • Date: 2005-03-16 02:15:29 UTC
  • Revision ID: james.westby@ubuntu.com-20050316021529-xrjlowsejs0cijig
Tags: upstream-0.3.2
Import upstream version 0.3.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/linalg/iterative/BiCGREVCOM.f.src
 
1
*
 
2
      SUBROUTINE <_c>BICGREVCOM( N, B, X, WORK, LDW, ITER, RESID, INFO,
 
3
     $                       NDX1, NDX2, SCLR1, SCLR2, IJOB)
 
4
*
 
5
*
 
6
*  -- Iterative template routine --
 
7
*     Univ. of Tennessee and Oak Ridge National Laboratory
 
8
*     October 1, 1993
 
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
*     Eijkhout, Pozo, Romine, and van der Vorst, SIAM Publications,
 
13
*     1993. (ftp netlib2.cs.utk.edu; cd linalg; get templates.ps).
 
14
*
 
15
*     .. Scalar Arguments ..
 
16
      INTEGER            N, LDW, ITER, INFO
 
17
      <rt=real,double precision,real,double precision>  RESID
 
18
      INTEGER            NDX1, NDX2
 
19
      <_t>   SCLR1, SCLR2
 
20
      INTEGER            IJOB
 
21
*     ..
 
22
*     .. Array Arguments ..
 
23
      <_t>   X( * ), B( * ), WORK( LDW,* )
 
24
*
 
25
*     ..
 
26
*
 
27
*  Purpose
 
28
*  =======
 
29
*
 
30
*  BiCG solves the linear system Ax = b using the
 
31
*  BiConjugate Gradient iterative method with preconditioning.
 
32
*
 
33
*  Arguments
 
34
*  =========
 
35
*
 
36
*  N       (input) INTEGER. 
 
37
*          On entry, the dimension of the matrix.
 
38
*          Unchanged on exit.
 
39
*
 
40
*  B       (input) DOUBLE PRECISION array, dimension N.
 
41
*          On entry, right hand side vector B.
 
42
*          Unchanged on exit.
 
43
*
 
44
*  X      (input/output) DOUBLE PRECISION array, dimension N.
 
45
*          On input, the initial guess; on exit, the iterated solution.
 
46
*
 
47
*  WORK    (workspace) DOUBLE PRECISION array, dimension (LDW,6).
 
48
*          Workspace for residual, direction vector, etc.
 
49
*          Note that Z and Q, and ZTLD and QTLD share workspace.
 
50
*
 
51
*  LDW     (input) INTEGER
 
52
*          The leading dimension of the array WORK. LDW .gt. = max(1,N).
 
53
*
 
54
*  ITER    (input/output) INTEGER
 
55
*          On input, the maximum iterations to be performed.
 
56
*          On output, actual number of iterations performed.
 
57
*
 
58
*  RESID   (input/output) DOUBLE PRECISION
 
59
*          On input, the allowable convergence measure for
 
60
*          norm( b - A*x ) / norm( b ).
 
61
*          On output, the final value of this measure.
 
62
*
 
63
*  INFO    (output) INTEGER
 
64
*
 
65
*          =  0: Successful exit. Iterated approximate solution returned.
 
66
*
 
67
*          .gt.   0: Convergence to tolerance not achieved. This will be
 
68
*                set to the number of iterations performed.
 
69
*
 
70
*          .ls.   0: Illegal input parameter, or breakdown occurred
 
71
*                during iteration.
 
72
*
 
73
*                Illegal parameter:
 
74
*
 
75
*                   -1: matrix dimension N .ls.  0
 
76
*                   -2: LDW .ls.  N
 
77
*                   -3: Maximum number of iterations ITER .ls. = 0.
 
78
*                   -5: Erroneous NDX1/NDX2 in INIT call.
 
79
*                   -6: Erroneous RLBL.
 
80
*
 
81
*                BREAKDOWN: If parameters RHO or OMEGA become smaller
 
82
*                   than some tolerance, the program will terminate.
 
83
*                   Here we check against tolerance BREAKTOL.
 
84
*
 
85
*                  -10: RHO .ls.  BREAKTOL: RHO and RTLD have become
 
86
*                                       orthogonal.
 
87
*
 
88
*                  BREAKTOL is set in func GETBREAK.
 
89
*
 
90
*  NDX1    (input/output) INTEGER. 
 
91
*  NDX2    On entry in INIT call contain indices required by interface
 
92
*          level for stopping test.
 
93
*          All other times, used as output, to indicate indices into
 
94
*          WORK[] for the MATVEC, PSOLVE done by the interface level.
 
95
*
 
96
*  SCLR1   (output) DOUBLE PRECISION.
 
97
*  SCLR2   Used to pass the scalars used in MATVEC. Scalars are reqd because
 
98
*          original routines use dgemv.
 
99
*
 
100
*  IJOB    (input/output) INTEGER. 
 
101
*          Used to communicate job code between the two levels.
 
102
*
 
103
*  BLAS CALLS:   DAXPY, DCOPY, DDOT, DNRM2,
 
104
*  ==============================================================
 
105
*
 
106
*     .. Parameters ..
 
107
      <rt>   ZERO, ONE
 
108
      PARAMETER        ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 
109
*     ..
 
110
*     .. Local Scalars ..
 
111
      INTEGER            R, RTLD, Z, ZTLD, P, PTLD, Q, QTLD, MAXIT,
 
112
     $                   NEED1, NEED2
 
113
      <rt>   TOL, BNRM2, RHOTOL, 
 
114
     $       <sdsd=s,d,s,d>GETBREAK,  
 
115
     $       <rc=s,d,sc,dz>NRM2
 
116
      <_t>  ALPHA, BETA, RHO, RHO1, <xdot=sdot,ddot,cdotc,zdotc>
 
117
*
 
118
*     indicates where to resume from. Only valid when IJOB = 2!
 
119
      INTEGER RLBL
 
120
*
 
121
*     saving all.
 
122
      SAVE
 
123
*
 
124
*     ..
 
125
*     .. External Routines ..
 
126
      EXTERNAL           <_c>AXPY, <_c>COPY, <xdot>, <rc>NRM2
 
127
*     ..
 
128
*     .. Executable Statements ..
 
129
*
 
130
*     Entry point, so test IJOB
 
131
      IF (IJOB .eq. 1) THEN
 
132
         GOTO 1
 
133
      ELSEIF (IJOB .eq. 2) THEN
 
134
*        here we do resumption handling
 
135
         IF (RLBL .eq. 2) GOTO 2
 
136
         IF (RLBL .eq. 3) GOTO 3
 
137
         IF (RLBL .eq. 4) GOTO 4
 
138
         IF (RLBL .eq. 5) GOTO 5
 
139
         IF (RLBL .eq. 6) GOTO 6
 
140
         IF (RLBL .eq. 7) GOTO 7
 
141
*        if neither of these, then error
 
142
         INFO = -6
 
143
         GOTO 20
 
144
      ENDIF
 
145
*
 
146
* init.
 
147
*****************
 
148
 1    CONTINUE
 
149
*****************
 
150
*
 
151
      INFO = 0
 
152
      MAXIT = ITER
 
153
      TOL   = RESID
 
154
*
 
155
*     Alias workspace columns.
 
156
*
 
157
      R    = 1
 
158
      RTLD = 2
 
159
      Z    = 3
 
160
      ZTLD = 4
 
161
      P    = 5
 
162
      PTLD = 6
 
163
      Q    = 3
 
164
      QTLD = 4
 
165
*
 
166
*     Check if caller will need indexing info.
 
167
*
 
168
      IF( NDX1.NE.-1 ) THEN
 
169
         IF( NDX1.EQ.1 ) THEN
 
170
            NEED1 = ((R - 1) * LDW) + 1
 
171
         ELSEIF( NDX1.EQ.2 ) THEN
 
172
            NEED1 = ((RTLD - 1) * LDW) + 1
 
173
         ELSEIF( NDX1.EQ.3 ) THEN
 
174
            NEED1 = ((Z - 1) * LDW) + 1
 
175
         ELSEIF( NDX1.EQ.4 ) THEN
 
176
            NEED1 = ((ZTLD - 1) * LDW) + 1
 
177
         ELSEIF( NDX1.EQ.5 ) THEN
 
178
            NEED1 = ((P - 1) * LDW) + 1
 
179
         ELSEIF( NDX1.EQ.6 ) THEN
 
180
            NEED1 = ((PTLD - 1) * LDW) + 1
 
181
         ELSEIF( NDX1.EQ.7 ) THEN
 
182
            NEED1 = ((Q - 1) * LDW) + 1
 
183
         ELSEIF( NDX1.EQ.8 ) THEN
 
184
            NEED1 = ((QTLD - 1) * LDW) + 1
 
185
         ELSE
 
186
*           report error
 
187
            INFO = -5
 
188
            GO TO 20
 
189
         ENDIF
 
190
      ELSE
 
191
         NEED1 = NDX1
 
192
      ENDIF
 
193
*
 
194
      IF( NDX2.NE.-1 ) THEN
 
195
         IF( NDX2.EQ.1 ) THEN
 
196
            NEED2 = ((R - 1) * LDW) + 1
 
197
         ELSEIF( NDX2.EQ.2 ) THEN
 
198
            NEED2 = ((RTLD - 1) * LDW) + 1
 
199
         ELSEIF( NDX2.EQ.3 ) THEN
 
200
            NEED2 = ((Z - 1) * LDW) + 1
 
201
         ELSEIF( NDX2.EQ.4 ) THEN
 
202
            NEED2 = ((ZTLD - 1) * LDW) + 1
 
203
         ELSEIF( NDX2.EQ.5 ) THEN
 
204
            NEED2 = ((P - 1) * LDW) + 1
 
205
         ELSEIF( NDX2.EQ.6 ) THEN
 
206
            NEED2 = ((PTLD - 1) * LDW) + 1
 
207
         ELSEIF( NDX2.EQ.7 ) THEN
 
208
            NEED2 = ((Q - 1) * LDW) + 1
 
209
         ELSEIF( NDX2.EQ.8 ) THEN
 
210
            NEED2 = ((QTLD - 1) * LDW) + 1
 
211
         ELSE
 
212
*           report error
 
213
            INFO = -5
 
214
            GO TO 20
 
215
         ENDIF
 
216
      ELSE
 
217
         NEED2 = NDX2
 
218
      ENDIF
 
219
*
 
220
*     Set breakdown parameters.
 
221
*
 
222
      RHOTOL = <sdsd>GETBREAK()
 
223
* 
 
224
*     Set initial residual.
 
225
*
 
226
      CALL <_c>COPY( N, B, 1, WORK(1,R), 1 )
 
227
      IF ( <rc>NRM2( N, X, 1 ).NE.ZERO ) THEN
 
228
*********CALL MATVEC( -ONE, X, ZERO, WORK(1,R) )
 
229
*        using WORK[RTLD] as temp
 
230
*********CALL <_c>COPY( N, X, 1, WORK(1,RTLD), 1 )
 
231
         SCLR1 = -ONE
 
232
         SCLR2 = ZERO
 
233
         NDX1 = ((RTLD - 1) * LDW) + 1
 
234
         NDX2 = ((R    - 1) * LDW) + 1
 
235
         RLBL = 2
 
236
         IJOB = 5
 
237
         RETURN
 
238
*****************
 
239
 2       CONTINUE
 
240
*****************
 
241
*
 
242
         IF ( <rc>NRM2( N, WORK(1,R), 1 ).LE.TOL ) GO TO 30
 
243
*
 
244
      ENDIF
 
245
      CALL <_c>COPY( N, WORK(1,R), 1, WORK(1,RTLD), 1 )
 
246
      BNRM2 = <rc>NRM2( N, B, 1 )
 
247
      IF ( BNRM2.EQ.ZERO ) BNRM2 = ONE
 
248
*
 
249
      ITER = 0
 
250
*
 
251
   10 CONTINUE
 
252
*
 
253
*     Perform BiConjugate Gradient iteration.
 
254
*
 
255
         ITER = ITER + 1
 
256
*
 
257
*        Compute direction vectors PK and PTLD.
 
258
*
 
259
*********CALL PSOLVE( WORK(1,Z), WORK(1,R) )
 
260
*
 
261
         NDX1 = ((Z - 1) * LDW) + 1
 
262
         NDX2 = ((R - 1) * LDW) + 1
 
263
         RLBL = 3
 
264
         IJOB = 3
 
265
         RETURN
 
266
*****************
 
267
 3       CONTINUE
 
268
*****************
 
269
*
 
270
*********CALL PSOLVETRANS( WORK(1,ZTLD), WORK(1,RTLD) )
 
271
*
 
272
         NDX1 = ((ZTLD - 1) * LDW) + 1
 
273
         NDX2 = ((RTLD - 1) * LDW) + 1
 
274
         RLBL = 4
 
275
         IJOB = 4
 
276
         RETURN
 
277
*****************
 
278
 4       CONTINUE
 
279
*****************
 
280
*
 
281
*
 
282
*         RHO = <xdot>( N, WORK(1,Z), 1, WORK(1,RTLD), 1 )
 
283
         RHO = <xdot>( N, WORK(1,RTLD), 1, WORK(1,Z), 1 )
 
284
         IF ( ABS( RHO ).LT.RHOTOL ) GO TO 25
 
285
*
 
286
         IF ( ITER.GT.1 ) THEN
 
287
            BETA = RHO / RHO1
 
288
            CALL <_c>AXPY( N, BETA, WORK(1,P), 1, WORK(1,Z), 1 )
 
289
*            CALL <_c>AXPY( N, BETA, WORK(1,PTLD), 1, WORK(1,ZTLD), 1 )
 
290
            CALL <_c>AXPY( N, <co= , ,conjg,conjg>(BETA), 
 
291
     $           WORK(1,PTLD), 1, WORK(1,ZTLD), 1 )
 
292
            CALL <_c>COPY( N, WORK(1,Z), 1, WORK(1,P), 1 )
 
293
            CALL <_c>COPY( N, WORK(1,ZTLD), 1, WORK(1,PTLD), 1 )
 
294
         ELSE
 
295
            CALL <_c>COPY( N, WORK(1,Z), 1, WORK(1,P), 1 )
 
296
            CALL <_c>COPY( N, WORK(1,ZTLD), 1, WORK(1,PTLD), 1 )
 
297
         ENDIF
 
298
*
 
299
*********CALL MATVEC( ONE, WORK(1,P), ZERO, WORK(1,Q) )
 
300
*
 
301
         SCLR1 = ONE
 
302
         SCLR2 = ZERO
 
303
         NDX1 = ((P - 1) * LDW) + 1
 
304
         NDX2 = ((Q - 1) * LDW) + 1
 
305
         RLBL = 5
 
306
         IJOB = 1
 
307
         RETURN
 
308
*****************
 
309
 5       CONTINUE
 
310
*****************
 
311
*
 
312
*********CALL MATVECTRANS( ONE, WORK(1,PTLD), ZERO, WORK(1,QTLD) )
 
313
*
 
314
         SCLR1 = ONE
 
315
         SCLR2 = ZERO
 
316
         NDX1 = ((PTLD - 1) * LDW) + 1
 
317
         NDX2 = ((QTLD - 1) * LDW) + 1
 
318
         RLBL = 6
 
319
         IJOB = 2
 
320
         RETURN
 
321
*****************
 
322
 6       CONTINUE
 
323
*****************
 
324
         ALPHA = RHO / <xdot>( N, WORK(1,PTLD), 1, WORK(1,Q), 1 )
 
325
*
 
326
*        Compute current solution vector x.
 
327
*
 
328
         CALL <_c>AXPY( N, ALPHA, WORK(1,P), 1, X, 1 )
 
329
*
 
330
*        Compute residual vector rk, find norm,
 
331
*        then check for tolerance.
 
332
*
 
333
         CALL <_c>AXPY( N, -ALPHA, WORK(1,Q), 1, WORK(1,R), 1 )
 
334
*
 
335
*********RESID = <rc>NRM2( N, WORK(1,R), 1 ) / BNRM2
 
336
*********IF ( RESID.LE.TOL  ) GO TO 30
 
337
*
 
338
         NDX1 = NEED1
 
339
         NDX2 = NEED2
 
340
*        Prepare for resumption & return
 
341
         RLBL = 7
 
342
         IJOB = 6
 
343
         RETURN
 
344
*
 
345
*****************
 
346
 7       CONTINUE
 
347
*****************
 
348
         IF( INFO.EQ.1 ) GO TO 30
 
349
*
 
350
         IF ( ITER.EQ.MAXIT ) THEN
 
351
            INFO = 1
 
352
            GO TO 20
 
353
         ENDIF
 
354
*
 
355
*         CALL <_c>AXPY( N, -ALPHA, WORK(1,QTLD), 1, WORK(1,RTLD), 1 )
 
356
         CALL <_c>AXPY( N, -<co>(ALPHA)
 
357
     $        , WORK(1,QTLD), 1, WORK(1,RTLD), 1 )
 
358
         RHO1 = RHO
 
359
*
 
360
      GO TO 10
 
361
*
 
362
   20 CONTINUE
 
363
*
 
364
*     Iteration fails.
 
365
*
 
366
      RLBL = -1
 
367
      IJOB = -1
 
368
      RETURN
 
369
*
 
370
   25 CONTINUE
 
371
*
 
372
*     Set breakdown flag.
 
373
*
 
374
      INFO = -10
 
375
      RLBL = -1
 
376
      IJOB = -1
 
377
      RETURN
 
378
*
 
379
   30 CONTINUE
 
380
*
 
381
*     Iteration successful; return.
 
382
*
 
383
      INFO = 0
 
384
      RLBL = -1
 
385
      IJOB = -1
 
386
      RETURN
 
387
*
 
388
*     End of BICGREVCOM
 
389
*
 
390
      END SUBROUTINE <_c>BICGREVCOM
 
391
 
 
392
 
 
393
 
 
394
 
 
395