CoCalc -- huti_tfqmr.F90

GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/fhutiter/src/huti_tfqmr.F90
³²⁰⁶ views
1
module huti_tfqmr
2
  use huti_aux
3
  implicit none
4
  !
5
  ! *
6
  ! *  Elmer, A Finite Element Software for Multiphysical Problems
7
  ! *
8
  ! *  Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland
9
  ! * 
10
  ! * This library is free software; you can redistribute it and/or
11
  ! * modify it under the terms of the GNU Lesser General Public
12
  ! * License as published by the Free Software Foundation; either
13
  ! * version 2.1 of the License, or (at your option) any later version.
14
  ! *
15
  ! * This library is distributed in the hope that it will be useful,
16
  ! * but WITHOUT ANY WARRANTY; without even the implied warranty of
17
  ! * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
18
  ! * Lesser General Public License for more details.
19
  ! * 
20
  ! * You should have received a copy of the GNU Lesser General Public
21
  ! * License along with this library (in file ../LGPL-2.1); if not, write 
22
  ! * to the Free Software Foundation, Inc., 51 Franklin Street, 
23
  ! * Fifth Floor, Boston, MA  02110-1301  USA
24
  ! *
25
  ! *****************************************************************************/
26

27
  !
28
  ! Subroutines to implement Transpose Free QMR iteration
29
  !
30

31
#include  "huti_fdefs.h" 
32

33
  !*************************************************************************
34
  !*************************************************************************
35
  !
36
  ! These subroutines are based on a paper by Roland W. Freund:
37
  ! "A Transpose-Free Quasi-Minimal Residual Algorithm for Non-Hermitian
38
  !  Linear Systems", 1993 (SIAM J. Sci. Comput, March 1993)
39
  !
40
  ! All matrix-vector operations are done externally, so we do not need
41
  ! to know about the matrix structure (sparse or dense). Memory allocation
42
  ! for the working arrays has also been done externally.
43

44
  !*************************************************************************
45
  ! Work array is used in the following order:
46
  ! work(:,1) = v
47
  ! work(:,2) = y
48
  ! work(:,3) = y new
49
  ! work(:,4) = r tilde (zero)
50
  ! work(:,5) = t1v (temporary for matrix-vector operations)
51
  ! work(:,6) = t2v (temporary for matrix-vector operations)
52
  ! work(:,7) = w
53
  ! work(:,8) = d
54
  ! work(:,9) = r
55
  ! work(:,10) = trv (temporary vector for residual computations)
56
  !
57
  !*************************************************************************
58
  ! Definitions to make the code more understandable and to make it look
59
  ! like the pseudo code
60
  !
61
#define  X  xvec 
62
#define  B  rhsvec 
63

64
#define  V  work(:,1) 
65
#define  V_ind  1 
66
#define  Y  work(:,2) 
67
#define  Y_ind  2 
68
#define  YNEW  work(:,3) 
69
#define  YNEW_ind  3 
70
#define  RTLD  work(:,4) 
71
#define  RTLD_ind  4 
72
#define  T1V  work(:,5) 
73
#define  T1V_ind  5 
74
#define  T2V  work(:,6) 
75
#define  T2V_ind  6 
76
#define  W  work(:,7) 
77
#define  W_ind  7 
78
#define  D  work(:,8) 
79
#define  D_ind  8 
80
#define  R  work(:,9) 
81
#define  R_ind  9 
82
#define  TRV  work(:,10) 
83
#define  TRV_ind  10 
84

85
  ! This is the magic ratio for upperb and tolerance used in upper bound
86
  ! convergence test
87
#define  UPPERB_TOL_RATIO  10.0 
88

89
contains
90

91
  !*************************************************************************
92
  !*************************************************************************
93
  ! Single precision version
94
  !*************************************************************************
95
  !*************************************************************************
96

97
  subroutine  huti_stfqmrsolv  ( ndim, wrkdim, xvec, rhsvec, ipar,&
98
       dpar, work, matvecsubr, pcondlsubr, pcondrsubr, &
99
       dotprodfun, normfun, stopcfun )
100

101
    use huti_interfaces
102
    implicit none
103

104
    procedure( mv_iface_s ), pointer :: matvecsubr
105
    procedure( pc_iface_s ), pointer :: pcondlsubr
106
    procedure( pc_iface_s ), pointer :: pcondrsubr
107
    procedure( dotp_iface_s ), pointer :: dotprodfun
108
    procedure( norm_iface_s ), pointer :: normfun
109
    procedure( stopc_iface_s ), pointer :: stopcfun
110

111
    ! Parameters
112

113
    integer :: ndim, wrkdim
114
    real, dimension(ndim) :: xvec, rhsvec
115
    integer, dimension(HUTI_IPAR_DFLTSIZE) :: ipar
116
    double precision, dimension(HUTI_DPAR_DFLTSIZE) :: dpar
117
    real, dimension(ndim,wrkdim) :: work
118

119
    ! Local variables
120

121
    real :: rho, oldrho, eta, tau, gamma, oldgamma, alpha
122
    real :: beta, c
123
    integer :: iter_count
124

125
    real :: residual, upperb, rhsnorm, precrhsnorm
126

127
    !
128
    ! End of variable declarations
129
    !*********************************************************************
130

131
    !*********************************************************************
132
    ! The actual TFQMR begins here (look the pseudo code in the
133
    ! "A Transpose-Free..."-paper, algorithm 5.1)
134
    !
135
    ! First the initialization part
136
    !
137

138
    iter_count = 1
139

140
    ! The following applies for all matrix operations in this solver
141

142
    HUTI_EXTOP_MATTYPE = HUTI_MAT_NOTTRPSED
143

144
    ! Norms of right-hand side vector are used in convergence tests
145

146
    if ( HUTI_STOPC .eq. HUTI_TRESID_SCALED_BYB .or. & 
147
         HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYB .or. &
148
         HUTI_STOPC .eq. HUTI_UPPERB_STOPC ) then
149
       rhsnorm = normfun( HUTI_NDIM, B, 1 )
150
    end if
151
    if ( HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYPRECB ) then
152
       call pcondlsubr( D, B, ipar )
153
       precrhsnorm = normfun( HUTI_NDIM, D, 1 )
154
    end if
155

156
    !
157
    ! Part 1A - 1C
158
    !
159

160
    ! Generate vector X if needed
161

162
    if ( HUTI_INITIALX .eq. HUTI_RANDOMX ) then
163
       call  huti_srandvec   ( X, ipar )
164
    else if ( HUTI_INITIALX .ne. HUTI_USERSUPPLIEDX ) then
165
       X = 1
166
    end if
167

168
    call pcondrsubr( D, X, ipar )
169
    call matvecsubr( D, R, ipar )
170
    D = B - R
171
    call pcondlsubr( R, D, ipar )
172

173
    Y = R; W = R
174
    call pcondrsubr( V, Y, ipar )
175
    call matvecsubr( V, D, ipar )
176
    call pcondlsubr( V, D, ipar )
177
    T2V = V
178

179
    D = 0
180
    tau = normfun( HUTI_NDIM, R, 1 )
181
    oldgamma = 0; gamma = 0; eta = 0
182

183
    RTLD = R
184
    oldrho = dotprodfun ( HUTI_NDIM, RTLD, 1, R, 1 )
185
    if ( oldrho .eq. 0 ) then
186
       HUTI_INFO = HUTI_TFQMR_RHO
187
       go to 1000
188
    end if
189

190
    !
191
    ! This is where the loop starts (that is we continue from here after
192
    ! the first iteration)
193
    !
194
    !
195
    ! Part 2A
196
    !
197

198
300 continue
199

200
    alpha = oldrho / dotprodfun( HUTI_NDIM, RTLD, 1, V, 1 )
201
    YNEW = Y - alpha * V
202

203
    !
204
    ! Part 2B
205
    !
206
    !
207
    ! This is the inner loop from 2n-1 to 2n
208
    !
209

210
    ! First the 2n-1 case
211

212
    ! Note: We have already MATRIX * Y in T2V
213

214
    W = W - alpha * T2V
215
    gamma = ( normfun( HUTI_NDIM, W, 1 )) / tau
216
    c = 1 / sqrt( 1 + gamma * gamma )
217
    tau = tau * gamma * c
218

219
    D = Y + ((oldgamma * oldgamma * eta) / alpha) * D
220
    eta = c * c * alpha
221
    X = X + eta * D
222

223
    oldgamma = gamma
224

225
    !
226
    ! Check the convergence against selected stopping criterion
227
    !
228

229
    select case (HUTI_STOPC)
230
    case (HUTI_TRUERESIDUAL)
231
       call pcondrsubr( TRV, X, ipar )
232
       call matvecsubr( TRV, R, ipar )
233
       TRV = R - B
234
       call pcondlsubr( R, TRV, ipar )
235
       residual = normfun( HUTI_NDIM, R, 1 )
236
    case (HUTI_TRESID_SCALED_BYB)
237
       call pcondrsubr( TRV, X, ipar )
238
       call matvecsubr( TRV, R, ipar )
239
       TRV = R - B
240
       call pcondlsubr( R, TRV, ipar )
241
       residual = normfun( HUTI_NDIM, R, 1 ) / rhsnorm
242
    case (HUTI_PSEUDORESIDUAL)
243
       call matvecsubr( X, R, ipar )
244
       R = R - B
245
       call pcondlsubr( TRV, R, ipar )
246
       residual = normfun( HUTI_NDIM, TRV, 1 )
247
    case (HUTI_PRESID_SCALED_BYB)
248
       call matvecsubr( X, R, ipar )
249
       R = R - B
250
       call pcondlsubr( TRV, R, ipar )
251
       residual = normfun( HUTI_NDIM, TRV, 1 ) / rhsnorm
252
    case (HUTI_PRESID_SCALED_BYPRECB)
253
       call matvecsubr( X, R, ipar )
254
       R = R - B
255
       call pcondlsubr( TRV, R, ipar )
256
       residual = normfun( HUTI_NDIM, TRV, 1 ) / precrhsnorm
257
    case (HUTI_XDIFF_NORM)
258
       R = eta * D
259
       residual = normfun( HUTI_NDIM, R, 1 )
260
    case (HUTI_UPPERB_STOPC)
261
       upperb = real( sqrt( 2.0 * iter_count ) * tau / rhsnorm)
262
       if ( ( upperb / HUTI_TOLERANCE ) .lt. UPPERB_TOL_RATIO ) then
263
          call pcondrsubr( TRV, X, ipar )
264
          call matvecsubr( TRV, R, ipar )
265
          TRV = R - B
266
          call pcondlsubr( R, TRV, ipar )
267
          residual = normfun( HUTI_NDIM, R, 1 ) / rhsnorm
268
       else
269
          residual = upperb
270
       end if
271
    case (HUTI_USUPPLIED_STOPC)
272
       residual = stopcfun( X, B, R, ipar, dpar )
273
    case default
274
       call pcondrsubr( TRV, X, ipar )
275
       call matvecsubr( TRV, R, ipar )
276
       TRV = R - B
277
       call pcondlsubr( R, TRV, ipar )
278
       residual = normfun( HUTI_NDIM, R, 1 )
279
    end select
280

281
    if ( residual .lt. HUTI_TOLERANCE ) then
282
       HUTI_INFO = HUTI_CONVERGENCE
283
       go to 1000
284
    end if
285

286
    IF( residual /= residual .OR. residual > HUTI_MAXTOLERANCE ) THEN
287
      HUTI_INFO = HUTI_DIVERGENCE
288
      GOTO 1000
289
    END IF
290

291
    !
292
    ! And then the 2n case
293
    !
294

295
    call pcondrsubr( T1V, YNEW, ipar )
296
    call matvecsubr( T1V, R, ipar )
297
    call pcondlsubr( T1V, R, ipar )
298

299
    W = W - alpha * T1V
300
    gamma = ( normfun( HUTI_NDIM, W, 1 )) / tau
301
    c = 1 / sqrt( 1 + gamma * gamma )
302
    tau = tau * gamma * c
303

304
    D = YNEW + ((oldgamma * oldgamma * eta) / alpha) * D
305
    eta = c * c * alpha
306
    X = X + eta * D
307

308
    oldgamma = gamma
309

310
    !
311
    ! Check the convergence against selected stopping criterion
312
    !
313

314
    select case (HUTI_STOPC)
315
    case (HUTI_TRUERESIDUAL)
316
       call pcondrsubr( TRV, X, ipar )
317
       call matvecsubr( TRV, R, ipar )
318
       TRV = R - B
319
       call pcondlsubr( R, TRV, ipar )
320
       residual = normfun( HUTI_NDIM, R, 1 )
321
    case (HUTI_TRESID_SCALED_BYB)
322
       call pcondrsubr( TRV, X, ipar )
323
       call matvecsubr( TRV, R, ipar )
324
       TRV = R - B
325
       call pcondlsubr( R, TRV, ipar )
326
       residual = normfun( HUTI_NDIM, R, 1 ) / rhsnorm
327
    case (HUTI_PSEUDORESIDUAL)
328
       call matvecsubr( X, R, ipar )
329
       R = R - B
330
       call pcondlsubr( TRV, R, ipar )
331
       residual = normfun( HUTI_NDIM, TRV, 1 )
332
    case (HUTI_PRESID_SCALED_BYB)
333
       call matvecsubr( X, R, ipar )
334
       R = R - B
335
       call pcondlsubr( TRV, R, ipar )
336
       residual = normfun( HUTI_NDIM, TRV, 1 ) / rhsnorm
337
    case (HUTI_PRESID_SCALED_BYPRECB)
338
       call matvecsubr( X, R, ipar )
339
       R = R - B
340
       call pcondlsubr( TRV, R, (/1/) )
341
       residual = normfun( HUTI_NDIM, TRV, 1 ) / precrhsnorm
342
    case (HUTI_XDIFF_NORM)
343
       R = eta * D
344
       residual = normfun( HUTI_NDIM, R, 1 )
345
    case (HUTI_UPPERB_STOPC)
346
       upperb = real( sqrt( 2.0 * iter_count ) * tau / rhsnorm)
347
       if ( ( upperb / HUTI_TOLERANCE ) .lt. UPPERB_TOL_RATIO ) then
348
          call pcondrsubr( TRV, X, ipar )
349
          call matvecsubr( TRV, R, ipar )
350
          TRV = R - B
351
          call pcondlsubr( R, TRV, ipar )
352
          residual = normfun( HUTI_NDIM, R, 1 ) / rhsnorm
353
       else
354
          residual = upperb
355
       end if
356
    case (HUTI_USUPPLIED_STOPC)
357
       residual = stopcfun( X, B, R, ipar, dpar )
358
    case default
359
       call pcondrsubr( TRV, X, ipar )
360
       call matvecsubr( TRV, R, ipar )
361
       TRV = R - B
362
       call pcondlsubr( R, TRV, ipar )
363
       residual = normfun( HUTI_NDIM, R, 1 )
364
    end select
365

366
    if ( residual .lt. HUTI_TOLERANCE ) then
367
       HUTI_INFO = HUTI_CONVERGENCE
368
       go to 1000
369
    end if
370

371
    IF( residual /= residual .OR. residual > HUTI_MAXTOLERANCE ) THEN
372
      HUTI_INFO = HUTI_DIVERGENCE
373
      GOTO 1000
374
    END IF
375

376
    !
377
    ! Produce debugging output if desired
378
    !
379

380
    if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
381
       if ( mod(iter_count, HUTI_DBUGLVL) .eq. 0 ) then
382
          if ( HUTI_STOPC .eq. HUTI_UPPERB_STOPC ) then
383
             write (*, '(I8, 2E17.7)') iter_count, residual, upperb
384
          else
385
             write (*, '(I8, E17.7)') iter_count, residual
386
          end if
387
       end if
388
    end if
389

390
    !
391
    ! Part 2C
392
    !
393

394
    rho = dotprodfun( HUTI_NDIM, RTLD, 1, W, 1 )
395
    beta = rho / oldrho
396
    YNEW = W + beta * YNEW
397
    call pcondrsubr( T2V, YNEW, ipar )
398
    call matvecsubr( T2V, R, ipar )
399
    call pcondlsubr( T2V, R, ipar )
400

401
    ! Note: we still have MATRIX * YNEW in T1V
402

403
    V = T2V + beta * T1V + beta * beta * V
404
    Y = YNEW
405

406
    oldrho = rho
407

408
    !
409
    ! Return next time back to the iteration loop (without initialization)
410
    !
411

412
    iter_count = iter_count + 1
413
    if ( iter_count .gt. HUTI_MAXIT ) then
414
       HUTI_INFO = HUTI_MAXITER
415
       go to 1000
416
    end if
417

418
    go to 300
419

420
    !
421
    ! This is where we exit last time (after enough iterations or breakdown)
422
    !
423

424
1000 continue
425

426
    ! Compute the unpreconditioned X
427

428
    call pcondrsubr( TRV, X, ipar )
429
    X = TRV
430

431
    if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
432
       if ( HUTI_STOPC .eq. HUTI_UPPERB_STOPC ) then
433
          write (*, '(I8, 2E17.7)') iter_count, residual, upperb
434
       else
435
          write (*, '(I8, E17.7)') iter_count, residual
436
       end if
437
    end if
438
    HUTI_ITERS = iter_count
439

440
    return
441

442
    ! End of execution
443
    !*********************************************************************
444

445
  end subroutine  huti_stfqmrsolv
446

447
  !*************************************************************************
448

449
  !*************************************************************************
450
  !*************************************************************************
451
  ! Double precision version
452
  !*************************************************************************
453
  !*************************************************************************
454

455
  subroutine  huti_dtfqmrsolv  ( ndim, wrkdim, xvec, rhsvec, ipar,&
456
       dpar, work, matvecsubr, pcondlsubr, pcondrsubr, &
457
       dotprodfun, normfun, stopcfun )
458

459
    use huti_interfaces
460
    implicit none
461

462
    procedure( mv_iface_d ), pointer :: matvecsubr
463
    procedure( pc_iface_d ), pointer :: pcondlsubr
464
    procedure( pc_iface_d ), pointer :: pcondrsubr
465
    procedure( dotp_iface_d ), pointer :: dotprodfun
466
    procedure( norm_iface_d ), pointer :: normfun
467
    procedure( stopc_iface_d ), pointer :: stopcfun
468

469
    ! Parameters
470

471
    integer :: ndim, wrkdim
472
    double precision, dimension(ndim) :: xvec, rhsvec
473
    integer, dimension(HUTI_IPAR_DFLTSIZE) :: ipar
474
    double precision, dimension(HUTI_DPAR_DFLTSIZE) :: dpar
475
    double precision, dimension(ndim,wrkdim) :: work
476

477
    ! Local variables
478

479
    double precision :: rho, oldrho, eta, tau, gamma, oldgamma, alpha
480
    double precision :: beta, c
481
    integer :: iter_count
482

483
    double precision :: residual, upperb, rhsnorm, precrhsnorm
484

485
    !
486
    ! End of variable declarations
487
    !*********************************************************************
488

489
    !*********************************************************************
490
    ! The actual TFQMR begins here (look the pseudo code in the
491
    ! "A Transpose-Free..."-paper, algorithm 5.1)
492
    !
493
    ! First the initialization part
494
    !
495

496
    iter_count = 1
497

498
    ! The following applies for all matrix operations in this solver
499

500
    HUTI_EXTOP_MATTYPE = HUTI_MAT_NOTTRPSED
501

502
    ! Norms of right-hand side vector are used in convergence tests
503

504
    if ( HUTI_STOPC .eq. HUTI_TRESID_SCALED_BYB .or. & 
505
         HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYB .or. &
506
         HUTI_STOPC .eq. HUTI_UPPERB_STOPC ) then
507
       rhsnorm = normfun( HUTI_NDIM, B, 1 )
508
    end if
509
    if ( HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYPRECB ) then
510
       call pcondlsubr( D, B, ipar )
511
       precrhsnorm = normfun( HUTI_NDIM, D, 1 )
512
    end if
513

514
    !
515
    ! Part 1A - 1C
516
    !
517

518
    ! Generate vector X if needed
519

520
    if ( HUTI_INITIALX .eq. HUTI_RANDOMX ) then
521
       call  huti_drandvec   ( X, ipar )
522
    else if ( HUTI_INITIALX .ne. HUTI_USERSUPPLIEDX ) then
523
       X = 1
524
    end if
525

526
    call pcondrsubr( D, X, ipar )
527
    call matvecsubr( D, R, ipar )
528
    D = B - R
529
    call pcondlsubr( R, D, ipar )
530

531
    Y = R; W = R
532
    call pcondrsubr( V, Y, ipar )
533
    call matvecsubr( V, D, ipar )
534
    call pcondlsubr( V, D, ipar )
535
    T2V = V
536

537
    D = 0
538
    tau = normfun( HUTI_NDIM, R, 1 )
539
    oldgamma = 0; gamma = 0; eta = 0
540

541
    RTLD = R
542
    oldrho = dotprodfun ( HUTI_NDIM, RTLD, 1, R, 1 )
543
    if ( oldrho .eq. 0 ) then
544
       HUTI_INFO = HUTI_TFQMR_RHO
545
       go to 1000
546
    end if
547

548
    !
549
    ! This is where the loop starts (that is we continue from here after
550
    ! the first iteration)
551
    !
552
    !
553
    ! Part 2A
554
    !
555

556
300 continue
557

558
    alpha = oldrho / dotprodfun( HUTI_NDIM, RTLD, 1, V, 1 )
559
    YNEW = Y - alpha * V
560

561
    !
562
    ! Part 2B
563
    !
564
    !
565
    ! This is the inner loop from 2n-1 to 2n
566
    !
567

568
    ! First the 2n-1 case
569

570
    ! Note: We have already MATRIX * Y in T2V
571

572
    W = W - alpha * T2V
573
    gamma = ( normfun( HUTI_NDIM, W, 1 )) / tau
574
    c = 1 / sqrt( 1 + gamma * gamma )
575
    tau = tau * gamma * c
576

577
    D = Y + ((oldgamma * oldgamma * eta) / alpha) * D
578
    eta = c * c * alpha
579
    X = X + eta * D
580

581
    oldgamma = gamma
582

583
    !
584
    ! Check the convergence against selected stopping criterion
585
    !
586

587
    select case (HUTI_STOPC)
588
    case (HUTI_TRUERESIDUAL)
589
       call pcondrsubr( TRV, X, ipar )
590
       call matvecsubr( TRV, R, ipar )
591
       TRV = R - B
592
       call pcondlsubr( R, TRV, ipar )
593
       residual = normfun( HUTI_NDIM, R, 1 )
594
    case (HUTI_TRESID_SCALED_BYB)
595
       call pcondrsubr( TRV, X, ipar )
596
       call matvecsubr( TRV, R, ipar )
597
       TRV = R - B
598
       call pcondlsubr( R, TRV, ipar )
599
       residual = normfun( HUTI_NDIM, R, 1 ) / rhsnorm
600
    case (HUTI_PSEUDORESIDUAL)
601
       call matvecsubr( X, R, ipar )
602
       R = R - B
603
       call pcondlsubr( TRV, R, ipar )
604
       residual = normfun( HUTI_NDIM, TRV, 1 )
605
    case (HUTI_PRESID_SCALED_BYB)
606
       call matvecsubr( X, R, ipar )
607
       R = R - B
608
       call pcondlsubr( TRV, R, ipar )
609
       residual = normfun( HUTI_NDIM, TRV, 1 ) / rhsnorm
610
    case (HUTI_PRESID_SCALED_BYPRECB)
611
       call matvecsubr( X, R, ipar )
612
       R = R - B
613
       call pcondlsubr( TRV, R, ipar )
614
       residual = normfun( HUTI_NDIM, TRV, 1 ) / precrhsnorm
615
    case (HUTI_XDIFF_NORM)
616
       R = eta * D
617
       residual = normfun( HUTI_NDIM, R, 1 )
618
    case (HUTI_UPPERB_STOPC)
619
       upperb = real( sqrt( 2.0 * iter_count ) * tau / rhsnorm)
620
       if ( ( upperb / HUTI_TOLERANCE ) .lt. UPPERB_TOL_RATIO ) then
621
          call pcondrsubr( TRV, X, ipar )
622
          call matvecsubr( TRV, R, ipar )
623
          TRV = R - B
624
          call pcondlsubr( R, TRV, ipar )
625
          residual = normfun( HUTI_NDIM, R, 1 ) / rhsnorm
626
       else
627
          residual = upperb
628
       end if
629
    case (HUTI_USUPPLIED_STOPC)
630
       residual = stopcfun( X, B, R, ipar, dpar )
631
    case default
632
       call pcondrsubr( TRV, X, ipar )
633
       call matvecsubr( TRV, R, ipar )
634
       TRV = R - B
635
       call pcondlsubr( R, TRV, ipar )
636
       residual = normfun( HUTI_NDIM, R, 1 )
637
    end select
638

639
    if ( residual .lt. HUTI_TOLERANCE ) then
640
       HUTI_INFO = HUTI_CONVERGENCE
641
       go to 1000
642
    end if
643

644
    IF( residual /= residual .OR. residual > HUTI_MAXTOLERANCE ) THEN
645
      HUTI_INFO = HUTI_DIVERGENCE
646
      GOTO 1000
647
    END IF
648

649
    !
650
    ! And then the 2n case
651
    !
652

653
    call pcondrsubr( T1V, YNEW, ipar )
654
    call matvecsubr( T1V, R, ipar )
655
    call pcondlsubr( T1V, R, ipar )
656

657
    W = W - alpha * T1V
658
    gamma = ( normfun( HUTI_NDIM, W, 1 )) / tau
659
    c = 1 / sqrt( 1 + gamma * gamma )
660
    tau = tau * gamma * c
661

662
    D = YNEW + ((oldgamma * oldgamma * eta) / alpha) * D
663
    eta = c * c * alpha
664
    X = X + eta * D
665

666
    oldgamma = gamma
667

668
    !
669
    ! Check the convergence against selected stopping criterion
670
    !
671

672
    select case (HUTI_STOPC)
673
    case (HUTI_TRUERESIDUAL)
674
       call pcondrsubr( TRV, X, ipar )
675
       call matvecsubr( TRV, R, ipar )
676
       TRV = R - B
677
       call pcondlsubr( R, TRV, ipar )
678
       residual = normfun( HUTI_NDIM, R, 1 )
679
    case (HUTI_TRESID_SCALED_BYB)
680
       call pcondrsubr( TRV, X, ipar )
681
       call matvecsubr( TRV, R, ipar )
682
       TRV = R - B
683
       call pcondlsubr( R, TRV, ipar )
684
       residual = normfun( HUTI_NDIM, R, 1 ) / rhsnorm
685
    case (HUTI_PSEUDORESIDUAL)
686
       call matvecsubr( X, R, ipar )
687
       R = R - B
688
       call pcondlsubr( TRV, R, ipar )
689
       residual = normfun( HUTI_NDIM, TRV, 1 )
690
    case (HUTI_PRESID_SCALED_BYB)
691
       call matvecsubr( X, R, ipar )
692
       R = R - B
693
       call pcondlsubr( TRV, R, ipar )
694
       residual = normfun( HUTI_NDIM, TRV, 1 ) / rhsnorm
695
    case (HUTI_PRESID_SCALED_BYPRECB)
696
       call matvecsubr( X, R, ipar )
697
       R = R - B
698
       call pcondlsubr( TRV, R, (/1/) )
699
       residual = normfun( HUTI_NDIM, TRV, 1 ) / precrhsnorm
700
    case (HUTI_XDIFF_NORM)
701
       R = eta * D
702
       residual = normfun( HUTI_NDIM, R, 1 )
703
    case (HUTI_UPPERB_STOPC)
704
       upperb = real( sqrt( 2.0 * iter_count ) * tau / rhsnorm)
705
       if ( ( upperb / HUTI_TOLERANCE ) .lt. UPPERB_TOL_RATIO ) then
706
          call pcondrsubr( TRV, X, ipar )
707
          call matvecsubr( TRV, R, ipar )
708
          TRV = R - B
709
          call pcondlsubr( R, TRV, ipar )
710
          residual = normfun( HUTI_NDIM, R, 1 ) / rhsnorm
711
       else
712
          residual = upperb
713
       end if
714
    case (HUTI_USUPPLIED_STOPC)
715
       residual = stopcfun( X, B, R, ipar, dpar )
716
    case default
717
       call pcondrsubr( TRV, X, ipar )
718
       call matvecsubr( TRV, R, ipar )
719
       TRV = R - B
720
       call pcondlsubr( R, TRV, ipar )
721
       residual = normfun( HUTI_NDIM, R, 1 )
722
    end select
723

724
    if ( residual .lt. HUTI_TOLERANCE ) then
725
       HUTI_INFO = HUTI_CONVERGENCE
726
       go to 1000
727
    end if
728

729
    IF( residual /= residual .OR. residual > HUTI_MAXTOLERANCE ) THEN
730
      HUTI_INFO = HUTI_DIVERGENCE
731
      GOTO 1000
732
    END IF
733

734
    !
735
    ! Produce debugging output if desired
736
    !
737

738
    if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
739
       if ( mod(iter_count, HUTI_DBUGLVL) .eq. 0 ) then
740
          if ( HUTI_STOPC .eq. HUTI_UPPERB_STOPC ) then
741
             write (*, '(I8, 2E17.7)') iter_count, residual, upperb
742
          else
743
             write (*, '(I8, E17.7)') iter_count, residual
744
          end if
745
       end if
746
    end if
747

748
    !
749
    ! Part 2C
750
    !
751

752
    rho = dotprodfun( HUTI_NDIM, RTLD, 1, W, 1 )
753
    beta = rho / oldrho
754
    YNEW = W + beta * YNEW
755
    call pcondrsubr( T2V, YNEW, ipar )
756
    call matvecsubr( T2V, R, ipar )
757
    call pcondlsubr( T2V, R, ipar )
758

759
    ! Note: we still have MATRIX * YNEW in T1V
760

761
    V = T2V + beta * T1V + beta * beta * V
762
    Y = YNEW
763

764
    oldrho = rho
765

766
    !
767
    ! Return next time back to the iteration loop (without initialization)
768
    !
769

770
    iter_count = iter_count + 1
771
    if ( iter_count .gt. HUTI_MAXIT ) then
772
       HUTI_INFO = HUTI_MAXITER
773
       go to 1000
774
    end if
775

776
    go to 300
777

778
    !
779
    ! This is where we exit last time (after enough iterations or breakdown)
780
    !
781

782
1000 continue
783

784
    ! Compute the unpreconditioned X
785

786
    call pcondrsubr( TRV, X, ipar )
787
    X = TRV
788

789
    if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
790
       if ( HUTI_STOPC .eq. HUTI_UPPERB_STOPC ) then
791
          write (*, '(I8, 2E17.7)') iter_count, residual, upperb
792
       else
793
          write (*, '(I8, E17.7)') iter_count, residual
794
       end if
795
    end if
796
    HUTI_ITERS = iter_count
797

798
    return
799

800
    ! End of execution
801
    !*********************************************************************
802

803
  end subroutine  huti_dtfqmrsolv
804

805
  !*************************************************************************
806

807
  !*************************************************************************
808
  !*************************************************************************
809
  ! Complex version
810
  !*************************************************************************
811
  !*************************************************************************
812

813
  subroutine  huti_ctfqmrsolv  ( ndim, wrkdim, xvec, rhsvec, ipar,&
814
       dpar, work, matvecsubr, pcondlsubr, pcondrsubr, &
815
       dotprodfun, normfun, stopcfun )
816

817
    use huti_interfaces
818
    implicit none
819

820
    procedure( mv_iface_c ), pointer :: matvecsubr
821
    procedure( pc_iface_c ), pointer :: pcondlsubr
822
    procedure( pc_iface_c ), pointer :: pcondrsubr
823
    procedure( dotp_iface_c ), pointer :: dotprodfun
824
    procedure( norm_iface_c ), pointer :: normfun
825
    procedure( stopc_iface_c ), pointer :: stopcfun
826

827
    ! Parameters
828

829
    integer :: ndim, wrkdim
830
    complex, dimension(ndim) :: xvec, rhsvec
831
    integer, dimension(HUTI_IPAR_DFLTSIZE) :: ipar
832
    double precision, dimension(HUTI_DPAR_DFLTSIZE) :: dpar
833
    complex, dimension(ndim,wrkdim) :: work
834

835
    ! Local variables
836

837
    complex :: rho, oldrho, eta, tau, gamma, oldgamma, alpha
838
    complex :: beta, c
839
    integer :: iter_count
840

841
    real :: residual, upperb, rhsnorm, precrhsnorm
842

843
    !
844
    ! End of variable declarations
845
    !*********************************************************************
846

847
    !*********************************************************************
848
    ! The actual TFQMR begins here (look the pseudo code in the
849
    ! "A Transpose-Free..."-paper, algorithm 5.1)
850
    !
851
    ! First the initialization part
852
    !
853

854
    iter_count = 1
855

856
    ! The following applies for all matrix operations in this solver
857

858
    HUTI_EXTOP_MATTYPE = HUTI_MAT_NOTTRPSED
859

860
    ! Norms of right-hand side vector are used in convergence tests
861

862
    if ( HUTI_STOPC .eq. HUTI_TRESID_SCALED_BYB .or. & 
863
         HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYB .or. &
864
         HUTI_STOPC .eq. HUTI_UPPERB_STOPC ) then
865
       rhsnorm = normfun( HUTI_NDIM, B, 1 )
866
    end if
867
    if ( HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYPRECB ) then
868
       call pcondlsubr( D, B, ipar )
869
       precrhsnorm = normfun( HUTI_NDIM, D, 1 )
870
    end if
871

872
    !
873
    ! Part 1A - 1C
874
    !
875

876
    ! Generate vector X if needed
877

878
    if ( HUTI_INITIALX .eq. HUTI_RANDOMX ) then
879
       call  huti_crandvec   ( X, ipar )
880
    else if ( HUTI_INITIALX .ne. HUTI_USERSUPPLIEDX ) then
881
       X = 1
882
    end if
883

884
    call pcondrsubr( D, X, ipar )
885
    call matvecsubr( D, R, ipar )
886
    D = B - R
887
    call pcondlsubr( R, D, ipar )
888

889
    Y = R; W = R
890
    call pcondrsubr( V, Y, ipar )
891
    call matvecsubr( V, D, ipar )
892
    call pcondlsubr( V, D, ipar )
893
    T2V = V
894

895
    D = 0
896
    tau = normfun( HUTI_NDIM, R, 1 )
897
    oldgamma = 0; gamma = 0; eta = 0
898

899
    RTLD = R
900
    oldrho = dotprodfun ( HUTI_NDIM, RTLD, 1, R, 1 )
901
    if ( oldrho .eq. 0 ) then
902
       HUTI_INFO = HUTI_TFQMR_RHO
903
       go to 1000
904
    end if
905

906
    !
907
    ! This is where the loop starts (that is we continue from here after
908
    ! the first iteration)
909
    !
910
    !
911
    ! Part 2A
912
    !
913

914
300 continue
915

916
    alpha = oldrho / dotprodfun( HUTI_NDIM, RTLD, 1, V, 1 )
917
    YNEW = Y - alpha * V
918

919
    !
920
    ! Part 2B
921
    !
922
    !
923
    ! This is the inner loop from 2n-1 to 2n
924
    !
925

926
    ! First the 2n-1 case
927

928
    ! Note: We have already MATRIX * Y in T2V
929

930
    W = W - alpha * T2V
931
    gamma = ( normfun( HUTI_NDIM, W, 1 )) / tau
932
    c = 1 / sqrt( 1 + gamma * gamma )
933
    tau = tau * gamma * c
934

935
    D = Y + ((oldgamma * oldgamma * eta) / alpha) * D
936
    eta = c * c * alpha
937
    X = X + eta * D
938

939
    oldgamma = gamma
940

941
    !
942
    ! Check the convergence against selected stopping criterion
943
    !
944

945
    select case (HUTI_STOPC)
946
    case (HUTI_TRUERESIDUAL)
947
       call pcondrsubr( TRV, X, ipar )
948
       call matvecsubr( TRV, R, ipar )
949
       TRV = R - B
950
       call pcondlsubr( R, TRV, ipar )
951
       residual = normfun( HUTI_NDIM, R, 1 )
952
    case (HUTI_TRESID_SCALED_BYB)
953
       call pcondrsubr( TRV, X, ipar )
954
       call matvecsubr( TRV, R, ipar )
955
       TRV = R - B
956
       call pcondlsubr( R, TRV, ipar )
957
       residual = normfun( HUTI_NDIM, R, 1 ) / rhsnorm
958
    case (HUTI_PSEUDORESIDUAL)
959
       call matvecsubr( X, R, ipar )
960
       R = R - B
961
       call pcondlsubr( TRV, R, ipar )
962
       residual = normfun( HUTI_NDIM, TRV, 1 )
963
    case (HUTI_PRESID_SCALED_BYB)
964
       call matvecsubr( X, R, ipar )
965
       R = R - B
966
       call pcondlsubr( TRV, R, ipar )
967
       residual = normfun( HUTI_NDIM, TRV, 1 ) / rhsnorm
968
    case (HUTI_PRESID_SCALED_BYPRECB)
969
       call matvecsubr( X, R, ipar )
970
       R = R - B
971
       call pcondlsubr( TRV, R, ipar )
972
       residual = normfun( HUTI_NDIM, TRV, 1 ) / precrhsnorm
973
    case (HUTI_XDIFF_NORM)
974
       R = eta * D
975
       residual = normfun( HUTI_NDIM, R, 1 )
976
    case (HUTI_UPPERB_STOPC)
977
       upperb = real( sqrt( 2.0 * iter_count ) * tau / rhsnorm)
978
       if ( ( upperb / HUTI_TOLERANCE ) .lt. UPPERB_TOL_RATIO ) then
979
          call pcondrsubr( TRV, X, ipar )
980
          call matvecsubr( TRV, R, ipar )
981
          TRV = R - B
982
          call pcondlsubr( R, TRV, ipar )
983
          residual = normfun( HUTI_NDIM, R, 1 ) / rhsnorm
984
       else
985
          residual = upperb
986
       end if
987
    case (HUTI_USUPPLIED_STOPC)
988
       residual = stopcfun( X, B, R, ipar, dpar )
989
    case default
990
       call pcondrsubr( TRV, X, ipar )
991
       call matvecsubr( TRV, R, ipar )
992
       TRV = R - B
993
       call pcondlsubr( R, TRV, ipar )
994
       residual = normfun( HUTI_NDIM, R, 1 )
995
    end select
996

997
    if ( residual .lt. HUTI_TOLERANCE ) then
998
       HUTI_INFO = HUTI_CONVERGENCE
999
       go to 1000
1000
    end if
1001

1002
    IF( residual /= residual .OR. residual > HUTI_MAXTOLERANCE ) THEN
1003
      HUTI_INFO = HUTI_DIVERGENCE
1004
      GOTO 1000
1005
    END IF
1006

1007
    !
1008
    ! And then the 2n case
1009
    !
1010

1011
    call pcondrsubr( T1V, YNEW, ipar )
1012
    call matvecsubr( T1V, R, ipar )
1013
    call pcondlsubr( T1V, R, ipar )
1014

1015
    W = W - alpha * T1V
1016
    gamma = ( normfun( HUTI_NDIM, W, 1 )) / tau
1017
    c = 1 / sqrt( 1 + gamma * gamma )
1018
    tau = tau * gamma * c
1019

1020
    D = YNEW + ((oldgamma * oldgamma * eta) / alpha) * D
1021
    eta = c * c * alpha
1022
    X = X + eta * D
1023

1024
    oldgamma = gamma
1025

1026
    !
1027
    ! Check the convergence against selected stopping criterion
1028
    !
1029

1030
    select case (HUTI_STOPC)
1031
    case (HUTI_TRUERESIDUAL)
1032
       call pcondrsubr( TRV, X, ipar )
1033
       call matvecsubr( TRV, R, ipar )
1034
       TRV = R - B
1035
       call pcondlsubr( R, TRV, ipar )
1036
       residual = normfun( HUTI_NDIM, R, 1 )
1037
    case (HUTI_TRESID_SCALED_BYB)
1038
       call pcondrsubr( TRV, X, ipar )
1039
       call matvecsubr( TRV, R, ipar )
1040
       TRV = R - B
1041
       call pcondlsubr( R, TRV, ipar )
1042
       residual = normfun( HUTI_NDIM, R, 1 ) / rhsnorm
1043
    case (HUTI_PSEUDORESIDUAL)
1044
       call matvecsubr( X, R, ipar )
1045
       R = R - B
1046
       call pcondlsubr( TRV, R, ipar )
1047
       residual = normfun( HUTI_NDIM, TRV, 1 )
1048
    case (HUTI_PRESID_SCALED_BYB)
1049
       call matvecsubr( X, R, ipar )
1050
       R = R - B
1051
       call pcondlsubr( TRV, R, ipar )
1052
       residual = normfun( HUTI_NDIM, TRV, 1 ) / rhsnorm
1053
    case (HUTI_PRESID_SCALED_BYPRECB)
1054
       call matvecsubr( X, R, ipar )
1055
       R = R - B
1056
       call pcondlsubr( TRV, R, (/1/) )
1057
       residual = normfun( HUTI_NDIM, TRV, 1 ) / precrhsnorm
1058
    case (HUTI_XDIFF_NORM)
1059
       R = eta * D
1060
       residual = normfun( HUTI_NDIM, R, 1 )
1061
    case (HUTI_UPPERB_STOPC)
1062
       upperb = real( sqrt( 2.0 * iter_count ) * tau / rhsnorm)
1063
       if ( ( upperb / HUTI_TOLERANCE ) .lt. UPPERB_TOL_RATIO ) then
1064
          call pcondrsubr( TRV, X, ipar )
1065
          call matvecsubr( TRV, R, ipar )
1066
          TRV = R - B
1067
          call pcondlsubr( R, TRV, ipar )
1068
          residual = normfun( HUTI_NDIM, R, 1 ) / rhsnorm
1069
       else
1070
          residual = upperb
1071
       end if
1072
    case (HUTI_USUPPLIED_STOPC)
1073
       residual = stopcfun( X, B, R, ipar, dpar )
1074
    case default
1075
       call pcondrsubr( TRV, X, ipar )
1076
       call matvecsubr( TRV, R, ipar )
1077
       TRV = R - B
1078
       call pcondlsubr( R, TRV, ipar )
1079
       residual = normfun( HUTI_NDIM, R, 1 )
1080
    end select
1081

1082
    if ( residual .lt. HUTI_TOLERANCE ) then
1083
       HUTI_INFO = HUTI_CONVERGENCE
1084
       go to 1000
1085
    end if
1086

1087
    IF( residual /= residual .OR. residual > HUTI_MAXTOLERANCE ) THEN
1088
      HUTI_INFO = HUTI_DIVERGENCE
1089
      GOTO 1000
1090
    END IF
1091

1092
    !
1093
    ! Produce debugging output if desired
1094
    !
1095

1096
    if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
1097
       if ( mod(iter_count, HUTI_DBUGLVL) .eq. 0 ) then
1098
          if ( HUTI_STOPC .eq. HUTI_UPPERB_STOPC ) then
1099
             write (*, '(I8, 2E17.7)') iter_count, residual, upperb
1100
          else
1101
             write (*, '(I8, E17.7)') iter_count, residual
1102
          end if
1103
       end if
1104
    end if
1105

1106
    !
1107
    ! Part 2C
1108
    !
1109

1110
    rho = dotprodfun( HUTI_NDIM, RTLD, 1, W, 1 )
1111
    beta = rho / oldrho
1112
    YNEW = W + beta * YNEW
1113
    call pcondrsubr( T2V, YNEW, ipar )
1114
    call matvecsubr( T2V, R, ipar )
1115
    call pcondlsubr( T2V, R, ipar )
1116

1117
    ! Note: we still have MATRIX * YNEW in T1V
1118

1119
    V = T2V + beta * T1V + beta * beta * V
1120
    Y = YNEW
1121

1122
    oldrho = rho
1123

1124
    !
1125
    ! Return next time back to the iteration loop (without initialization)
1126
    !
1127

1128
    iter_count = iter_count + 1
1129
    if ( iter_count .gt. HUTI_MAXIT ) then
1130
       HUTI_INFO = HUTI_MAXITER
1131
       go to 1000
1132
    end if
1133

1134
    go to 300
1135

1136
    !
1137
    ! This is where we exit last time (after enough iterations or breakdown)
1138
    !
1139

1140
1000 continue
1141

1142
    ! Compute the unpreconditioned X
1143

1144
    call pcondrsubr( TRV, X, ipar )
1145
    X = TRV
1146

1147
    if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
1148
       if ( HUTI_STOPC .eq. HUTI_UPPERB_STOPC ) then
1149
          write (*, '(I8, 2E17.7)') iter_count, residual, upperb
1150
       else
1151
          write (*, '(I8, E17.7)') iter_count, residual
1152
       end if
1153
    end if
1154
    HUTI_ITERS = iter_count
1155

1156
    return
1157

1158
    ! End of execution
1159
    !*********************************************************************
1160

1161
  end subroutine  huti_ctfqmrsolv
1162

1163
  !*************************************************************************
1164

1165
  !*************************************************************************
1166
  !*************************************************************************
1167
  ! Double complex version
1168
  !*************************************************************************
1169
  !*************************************************************************
1170

1171
  subroutine  huti_ztfqmrsolv  ( ndim, wrkdim, xvec, rhsvec, ipar,&
1172
       dpar, work, matvecsubr, pcondlsubr, pcondrsubr, &
1173
       dotprodfun, normfun, stopcfun )
1174

1175
    use huti_interfaces
1176
    implicit none
1177

1178
    procedure( mv_iface_z ), pointer :: matvecsubr
1179
    procedure( pc_iface_z ), pointer :: pcondlsubr
1180
    procedure( pc_iface_z ), pointer :: pcondrsubr
1181
    procedure( dotp_iface_z ), pointer :: dotprodfun
1182
    procedure( norm_iface_z ), pointer :: normfun
1183
    procedure( stopc_iface_z ), pointer :: stopcfun
1184

1185
    ! Parameters
1186

1187
    integer :: ndim, wrkdim
1188
    double complex, dimension(ndim) :: xvec, rhsvec
1189
    integer, dimension(HUTI_IPAR_DFLTSIZE) :: ipar
1190
    double precision, dimension(HUTI_DPAR_DFLTSIZE) :: dpar
1191
    double complex, dimension(ndim,wrkdim) :: work
1192

1193
    ! Local variables
1194

1195
    double complex :: rho, oldrho, eta, tau, gamma, oldgamma, alpha
1196
    double complex :: beta, c
1197
    integer :: iter_count
1198

1199
    double precision :: residual, upperb, rhsnorm, precrhsnorm
1200

1201
    !
1202
    ! End of variable declarations
1203
    !*********************************************************************
1204

1205
    !*********************************************************************
1206
    ! The actual TFQMR begins here (look the pseudo code in the
1207
    ! "A Transpose-Free..."-paper, algorithm 5.1)
1208
    !
1209
    ! First the initialization part
1210
    !
1211

1212
    iter_count = 1
1213

1214
    ! The following applies for all matrix operations in this solver
1215

1216
    HUTI_EXTOP_MATTYPE = HUTI_MAT_NOTTRPSED
1217

1218
    ! Norms of right-hand side vector are used in convergence tests
1219

1220
    if ( HUTI_STOPC .eq. HUTI_TRESID_SCALED_BYB .or. & 
1221
         HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYB .or. &
1222
         HUTI_STOPC .eq. HUTI_UPPERB_STOPC ) then
1223
       rhsnorm = normfun( HUTI_NDIM, B, 1 )
1224
    end if
1225
    if ( HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYPRECB ) then
1226
       call pcondlsubr( D, B, ipar )
1227
       precrhsnorm = normfun( HUTI_NDIM, D, 1 )
1228
    end if
1229

1230
    !
1231
    ! Part 1A - 1C
1232
    !
1233

1234
    ! Generate vector X if needed
1235

1236
    if ( HUTI_INITIALX .eq. HUTI_RANDOMX ) then
1237
       call  huti_zrandvec   ( X, ipar )
1238
    else if ( HUTI_INITIALX .ne. HUTI_USERSUPPLIEDX ) then
1239
       X = 1
1240
    end if
1241

1242
    call pcondrsubr( D, X, ipar )
1243
    call matvecsubr( D, R, ipar )
1244
    D = B - R
1245
    call pcondlsubr( R, D, ipar )
1246

1247
    Y = R; W = R
1248
    call pcondrsubr( V, Y, ipar )
1249
    call matvecsubr( V, D, ipar )
1250
    call pcondlsubr( V, D, ipar )
1251
    T2V = V
1252

1253
    D = 0
1254
    tau = normfun( HUTI_NDIM, R, 1 )
1255
    oldgamma = 0; gamma = 0; eta = 0
1256

1257
    RTLD = R
1258
    oldrho = dotprodfun ( HUTI_NDIM, RTLD, 1, R, 1 )
1259
    if ( oldrho .eq. 0 ) then
1260
       HUTI_INFO = HUTI_TFQMR_RHO
1261
       go to 1000
1262
    end if
1263

1264
    !
1265
    ! This is where the loop starts (that is we continue from here after
1266
    ! the first iteration)
1267
    !
1268
    !
1269
    ! Part 2A
1270
    !
1271

1272
300 continue
1273

1274
    alpha = oldrho / dotprodfun( HUTI_NDIM, RTLD, 1, V, 1 )
1275
    YNEW = Y - alpha * V
1276

1277
    !
1278
    ! Part 2B
1279
    !
1280
    !
1281
    ! This is the inner loop from 2n-1 to 2n
1282
    !
1283

1284
    ! First the 2n-1 case
1285

1286
    ! Note: We have already MATRIX * Y in T2V
1287

1288
    W = W - alpha * T2V
1289
    gamma = ( normfun( HUTI_NDIM, W, 1 )) / tau
1290
    c = 1 / sqrt( 1 + gamma * gamma )
1291
    tau = tau * gamma * c
1292

1293
    D = Y + ((oldgamma * oldgamma * eta) / alpha) * D
1294
    eta = c * c * alpha
1295
    X = X + eta * D
1296

1297
    oldgamma = gamma
1298

1299
    !
1300
    ! Check the convergence against selected stopping criterion
1301
    !
1302

1303
    select case (HUTI_STOPC)
1304
    case (HUTI_TRUERESIDUAL)
1305
       call pcondrsubr( TRV, X, ipar )
1306
       call matvecsubr( TRV, R, ipar )
1307
       TRV = R - B
1308
       call pcondlsubr( R, TRV, ipar )
1309
       residual = normfun( HUTI_NDIM, R, 1 )
1310
    case (HUTI_TRESID_SCALED_BYB)
1311
       call pcondrsubr( TRV, X, ipar )
1312
       call matvecsubr( TRV, R, ipar )
1313
       TRV = R - B
1314
       call pcondlsubr( R, TRV, ipar )
1315
       residual = normfun( HUTI_NDIM, R, 1 ) / rhsnorm
1316
    case (HUTI_PSEUDORESIDUAL)
1317
       call matvecsubr( X, R, ipar )
1318
       R = R - B
1319
       call pcondlsubr( TRV, R, ipar )
1320
       residual = normfun( HUTI_NDIM, TRV, 1 )
1321
    case (HUTI_PRESID_SCALED_BYB)
1322
       call matvecsubr( X, R, ipar )
1323
       R = R - B
1324
       call pcondlsubr( TRV, R, ipar )
1325
       residual = normfun( HUTI_NDIM, TRV, 1 ) / rhsnorm
1326
    case (HUTI_PRESID_SCALED_BYPRECB)
1327
       call matvecsubr( X, R, ipar )
1328
       R = R - B
1329
       call pcondlsubr( TRV, R, ipar )
1330
       residual = normfun( HUTI_NDIM, TRV, 1 ) / precrhsnorm
1331
    case (HUTI_XDIFF_NORM)
1332
       R = eta * D
1333
       residual = normfun( HUTI_NDIM, R, 1 )
1334
    case (HUTI_UPPERB_STOPC)
1335
       upperb = real( sqrt( 2.0 * iter_count ) * tau / rhsnorm)
1336
       if ( ( upperb / HUTI_TOLERANCE ) .lt. UPPERB_TOL_RATIO ) then
1337
          call pcondrsubr( TRV, X, ipar )
1338
          call matvecsubr( TRV, R, ipar )
1339
          TRV = R - B
1340
          call pcondlsubr( R, TRV, ipar )
1341
          residual = normfun( HUTI_NDIM, R, 1 ) / rhsnorm
1342
       else
1343
          residual = upperb
1344
       end if
1345
    case (HUTI_USUPPLIED_STOPC)
1346
       residual = stopcfun( X, B, R, ipar, dpar )
1347
    case default
1348
       call pcondrsubr( TRV, X, ipar )
1349
       call matvecsubr( TRV, R, ipar )
1350
       TRV = R - B
1351
       call pcondlsubr( R, TRV, ipar )
1352
       residual = normfun( HUTI_NDIM, R, 1 )
1353
    end select
1354

1355
    if ( residual .lt. HUTI_TOLERANCE ) then
1356
       HUTI_INFO = HUTI_CONVERGENCE
1357
       go to 1000
1358
    end if
1359

1360
    IF( residual /= residual .OR. residual > HUTI_MAXTOLERANCE ) THEN
1361
      HUTI_INFO = HUTI_DIVERGENCE
1362
      GOTO 1000
1363
    END IF
1364

1365
    !
1366
    ! And then the 2n case
1367
    !
1368

1369
    call pcondrsubr( T1V, YNEW, ipar )
1370
    call matvecsubr( T1V, R, ipar )
1371
    call pcondlsubr( T1V, R, ipar )
1372

1373
    W = W - alpha * T1V
1374
    gamma = ( normfun( HUTI_NDIM, W, 1 )) / tau
1375
    c = 1 / sqrt( 1 + gamma * gamma )
1376
    tau = tau * gamma * c
1377

1378
    D = YNEW + ((oldgamma * oldgamma * eta) / alpha) * D
1379
    eta = c * c * alpha
1380
    X = X + eta * D
1381

1382
    oldgamma = gamma
1383

1384
    !
1385
    ! Check the convergence against selected stopping criterion
1386
    !
1387

1388
    select case (HUTI_STOPC)
1389
    case (HUTI_TRUERESIDUAL)
1390
       call pcondrsubr( TRV, X, ipar )
1391
       call matvecsubr( TRV, R, ipar )
1392
       TRV = R - B
1393
       call pcondlsubr( R, TRV, ipar )
1394
       residual = normfun( HUTI_NDIM, R, 1 )
1395
    case (HUTI_TRESID_SCALED_BYB)
1396
       call pcondrsubr( TRV, X, ipar )
1397
       call matvecsubr( TRV, R, ipar )
1398
       TRV = R - B
1399
       call pcondlsubr( R, TRV, ipar )
1400
       residual = normfun( HUTI_NDIM, R, 1 ) / rhsnorm
1401
    case (HUTI_PSEUDORESIDUAL)
1402
       call matvecsubr( X, R, ipar )
1403
       R = R - B
1404
       call pcondlsubr( TRV, R, ipar )
1405
       residual = normfun( HUTI_NDIM, TRV, 1 )
1406
    case (HUTI_PRESID_SCALED_BYB)
1407
       call matvecsubr( X, R, ipar )
1408
       R = R - B
1409
       call pcondlsubr( TRV, R, ipar )
1410
       residual = normfun( HUTI_NDIM, TRV, 1 ) / rhsnorm
1411
    case (HUTI_PRESID_SCALED_BYPRECB)
1412
       call matvecsubr( X, R, ipar )
1413
       R = R - B
1414
       call pcondlsubr( TRV, R, (/1/) )
1415
       residual = normfun( HUTI_NDIM, TRV, 1 ) / precrhsnorm
1416
    case (HUTI_XDIFF_NORM)
1417
       R = eta * D
1418
       residual = normfun( HUTI_NDIM, R, 1 )
1419
    case (HUTI_UPPERB_STOPC)
1420
       upperb = real( sqrt( 2.0 * iter_count ) * tau / rhsnorm)
1421
       if ( ( upperb / HUTI_TOLERANCE ) .lt. UPPERB_TOL_RATIO ) then
1422
          call pcondrsubr( TRV, X, ipar )
1423
          call matvecsubr( TRV, R, ipar )
1424
          TRV = R - B
1425
          call pcondlsubr( R, TRV, ipar )
1426
          residual = normfun( HUTI_NDIM, R, 1 ) / rhsnorm
1427
       else
1428
          residual = upperb
1429
       end if
1430
    case (HUTI_USUPPLIED_STOPC)
1431
       residual = stopcfun( X, B, R, ipar, dpar )
1432
    case default
1433
       call pcondrsubr( TRV, X, ipar )
1434
       call matvecsubr( TRV, R, ipar )
1435
       TRV = R - B
1436
       call pcondlsubr( R, TRV, ipar )
1437
       residual = normfun( HUTI_NDIM, R, 1 )
1438
    end select
1439

1440
    if ( residual .lt. HUTI_TOLERANCE ) then
1441
       HUTI_INFO = HUTI_CONVERGENCE
1442
       go to 1000
1443
    end if
1444

1445
    IF( residual /= residual .OR. residual > HUTI_MAXTOLERANCE ) THEN
1446
      HUTI_INFO = HUTI_DIVERGENCE
1447
      GOTO 1000
1448
    END IF
1449

1450
    !
1451
    ! Produce debugging output if desired
1452
    !
1453

1454
    if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
1455
       if ( mod(iter_count, HUTI_DBUGLVL) .eq. 0 ) then
1456
          if ( HUTI_STOPC .eq. HUTI_UPPERB_STOPC ) then
1457
             write (*, '(I8, 2E17.7)') iter_count, residual, upperb
1458
          else
1459
             write (*, '(I8, E17.7)') iter_count, residual
1460
          end if
1461
       end if
1462
    end if
1463

1464
    !
1465
    ! Part 2C
1466
    !
1467

1468
    rho = dotprodfun( HUTI_NDIM, RTLD, 1, W, 1 )
1469
    beta = rho / oldrho
1470
    YNEW = W + beta * YNEW
1471
    call pcondrsubr( T2V, YNEW, ipar )
1472
    call matvecsubr( T2V, R, ipar )
1473
    call pcondlsubr( T2V, R, ipar )
1474

1475
    ! Note: we still have MATRIX * YNEW in T1V
1476

1477
    V = T2V + beta * T1V + beta * beta * V
1478
    Y = YNEW
1479

1480
    oldrho = rho
1481

1482
    !
1483
    ! Return next time back to the iteration loop (without initialization)
1484
    !
1485

1486
    iter_count = iter_count + 1
1487
    if ( iter_count .gt. HUTI_MAXIT ) then
1488
       HUTI_INFO = HUTI_MAXITER
1489
       go to 1000
1490
    end if
1491

1492
    go to 300
1493

1494
    !
1495
    ! This is where we exit last time (after enough iterations or breakdown)
1496
    !
1497

1498
1000 continue
1499

1500
    ! Compute the unpreconditioned X
1501

1502
    call pcondrsubr( TRV, X, ipar )
1503
    X = TRV
1504

1505
    if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
1506
       if ( HUTI_STOPC .eq. HUTI_UPPERB_STOPC ) then
1507
          write (*, '(I8, 2E17.7)') iter_count, residual, upperb
1508
       else
1509
          write (*, '(I8, E17.7)') iter_count, residual
1510
       end if
1511
    end if
1512
    HUTI_ITERS = iter_count
1513

1514
    return
1515

1516
    ! End of execution
1517
    !*********************************************************************
1518

1519
  end subroutine  huti_ztfqmrsolv
1520

1521
  !*************************************************************************
1522

1523
end module huti_tfqmr
1524

1525
Product

Resources

Company
