CoCalc -- huti_bicgstab

GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/fhutiter/src/huti_bicgstab_2.F90
³²⁰⁶ views
1
module huti_bicgstab_2
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
! Subroutine to implement BiConjugate Gradient Stabilised (2) iteration
29

30
#include  "huti_fdefs.h" 
31

32
!*************************************************************************
33
!*************************************************************************
34
!
35
! This subroutine is based on a paper by Henk A. Van der Vorst:
36
! "Parallel Iterative Solution Methods for Linear Systems arising from
37
!  Discretized PDE's". This is the Bi-CGSTAB(2) version.
38
!
39
! All matrix-vector operations are done externally, so we do not need
40
! to know about the matrix structure (sparse or dense). Memory allocation
41
! for the working arrays has also been done externally.
42

43
!*************************************************************************
44
! Work array is used in the following order:
45
! work(:,1) = r tilde (zero)
46
! work(:,2) = u
47
! work(:,3) = t1v
48
! work(:,4) = v
49
! work(:,5) = s
50
! work(:,6) = w
51
! work(:,7) = t
52
! work(:,8) = r
53
!
54
!*************************************************************************
55
! Definitions to make the code more understandable and to make it look
56
! like the pseudo code
57
!
58

59
#define  X  xvec 
60
#define  B  rhsvec 
61

62
#define  RTLD  work(:,1) 
63
#define  RTLD_ind  1 
64
#define  U  work(:,2) 
65
#define  U_ind  2 
66
#define  T1V  work(:,3) 
67
#define  T1V_ind  3 
68
#define  V  work(:,4) 
69
#define  V_ind  4 
70
#define  S  work(:,5) 
71
#define  S_ind  5 
72
#define  W  work(:,6) 
73
#define  W_ind  6 
74
#define  T  work(:,7) 
75
#define  T_ind  7 
76
#define  R  work(:,8) 
77
#define  R_ind  8 
78

79
!*************************************************************************
80

81
contains
82

83
  !*************************************************************************
84
  !*************************************************************************
85
  ! Single precision version
86
  !*************************************************************************
87
  !*************************************************************************
88

89
  subroutine  huti_sbicgstab_2solv  ( ndim, wrkdim, xvec, rhsvec, &
90
       ipar, dpar, work, matvecsubr, pcondlsubr, &
91
       pcondrsubr, dotprodfun, normfun, stopcfun )
92

93
    use huti_interfaces
94
    implicit none
95

96
    procedure( mv_iface_s ), pointer :: matvecsubr
97
    procedure( pc_iface_s ), pointer :: pcondlsubr 
98
    procedure( pc_iface_s ), pointer :: pcondrsubr 
99
    procedure( dotp_iface_s ), pointer :: dotprodfun 
100
    procedure( norm_iface_s ), pointer :: normfun 
101
    procedure( stopc_iface_s ), pointer :: stopcfun 
102

103
    ! Parameters
104

105
    integer :: ndim, wrkdim
106
    real, dimension(ndim) :: xvec, rhsvec
107
    integer, dimension(HUTI_IPAR_DFLTSIZE) :: ipar
108
    double precision, dimension(HUTI_DPAR_DFLTSIZE) :: dpar
109
    real, dimension(ndim,wrkdim) :: work
110

111
    ! Local variables
112

113
    real :: rho, oldrho, alpha, beta, omega1, omega2
114
    real :: tau, delta, myy
115
    integer :: iter_count
116

117
    real :: residual, rhsnorm, precrhsnorm
118

119
    !
120
    ! End of variable declarations
121
    !*********************************************************************
122

123
    !*********************************************************************
124
    ! The actual BiCGSTAB begins here (look the pseudo code in the
125
    ! "Templates..."-book, page 27)
126
    !
127
    ! First the initialization part
128
    !
129

130
    iter_count = 1
131

132
    ! The following applies for all matrix operations in this solver
133

134
    HUTI_EXTOP_MATTYPE = HUTI_MAT_NOTTRPSED
135

136
    ! Norms of right-hand side vector are used in convergence tests
137

138
    if ( HUTI_STOPC .eq. HUTI_TRESID_SCALED_BYB .or. & 
139
         HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYB ) then
140
       rhsnorm = normfun( HUTI_NDIM, B, 1 )
141
    end if
142
    if ( HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYPRECB ) then
143
       call pcondlsubr( T1V, B, ipar )
144
       precrhsnorm = normfun( HUTI_NDIM, T1V, 1 )
145
    end if
146

147
    ! Generate vector X if needed
148

149
    if ( HUTI_INITIALX .eq. HUTI_RANDOMX ) then
150
       call  huti_srandvec   ( X, ipar )
151
    else if ( HUTI_INITIALX .ne. HUTI_USERSUPPLIEDX ) then
152
       X = 1
153
    end if
154

155
    call pcondrsubr( U, X, ipar )
156
    call matvecsubr( U, R, ipar )
157
    U = B - R
158
    call pcondlsubr( R, U, ipar )
159
    RTLD = R
160
    U = 0
161
    oldrho = 1; omega2 = 1; alpha = 0
162

163
    !
164
    ! This is where the loop starts (that is we continue from here after
165
    ! the first iteration)
166
    !
167

168
300 continue
169

170
    oldrho = -omega2 * oldrho
171

172
    !
173
    ! This is the even BiCG step
174
    !
175

176
    rho = dotprodfun( HUTI_NDIM, RTLD, 1, R, 1 )
177
    if ( rho .eq. 0 ) then
178
       HUTI_INFO = HUTI_BICGSTAB_2_RHO
179
       go to 1000
180
    end if
181

182
    beta = ( rho * alpha ) / oldrho
183
    oldrho = rho
184
    U = R - beta * U
185

186
    call pcondrsubr( V, U, ipar )
187
    call matvecsubr( V, T1V, ipar )
188
    call pcondlsubr( V, T1V, ipar )
189

190
    alpha = oldrho / dotprodfun( HUTI_NDIM, RTLD, 1, V, 1 )
191
    R = R - alpha * V
192

193
    call pcondrsubr( S, R, ipar )
194
    call matvecsubr( S, T1V, ipar )
195
    call pcondlsubr( S, T1V, ipar )
196

197
    X = X + alpha * U
198

199
    !
200
    ! This is the odd BiCG step
201
    !
202

203
    rho = dotprodfun( HUTI_NDIM, RTLD, 1, S, 1 )
204
    if ( rho .eq. 0 ) then
205
       HUTI_INFO = HUTI_BICGSTAB_2_RHO
206
       go to 1000
207
    end if
208

209
    beta = ( rho * alpha ) / oldrho
210
    oldrho = rho
211
    V = S - beta * V
212

213
    call pcondrsubr( W, V, ipar )
214
    call matvecsubr( W, T1V, ipar )
215
    call pcondlsubr( W, T1V, ipar )
216

217
    alpha = oldrho / dotprodfun( HUTI_NDIM, RTLD, 1, W, 1 )
218
    U = R - beta * U
219
    R = R - alpha * V
220
    S = S - alpha * W
221

222
    call pcondrsubr( T, S, ipar )
223
    call matvecsubr( T, T1V, ipar )
224
    call pcondlsubr( T, T1V, ipar )
225

226
    !
227
    ! This is the GCR(2) part
228
    !
229

230
    omega1 = dotprodfun( HUTI_NDIM, R, 1, S, 1 )
231
    myy = dotprodfun( HUTI_NDIM, S, 1, S, 1 )
232
    delta = dotprodfun( HUTI_NDIM, S, 1, T, 1 )
233
    tau = dotprodfun( HUTI_NDIM, T, 1, T, 1 )
234
    omega2 = dotprodfun( HUTI_NDIM, R, 1, T, 1 )
235

236
    tau = tau - ( delta * delta ) / myy
237
    omega2 = ( omega2 - ( delta * omega1 ) / myy ) / tau
238
    omega1 = ( omega1 - delta * omega2 ) / myy
239

240
    X = X + omega1 * R + omega2 * S + alpha * U
241
    R = R - omega1 * S - omega2 * T
242

243
    !
244
    ! Check the convergence against selected stopping criterion
245
    !
246

247
    select case (HUTI_STOPC)
248
    case (HUTI_TRUERESIDUAL)
249
       call pcondrsubr( S, X, ipar )
250
       call matvecsubr( S, T1V, ipar )
251
       T1V = T1V - B
252
       call pcondlsubr( S, T1V, ipar )
253
       residual = normfun( HUTI_NDIM, S, 1 )
254
    case (HUTI_TRESID_SCALED_BYB)
255
       call pcondrsubr( S, X, ipar )
256
       call matvecsubr( S, T1V, ipar )
257
       T1V = T1V - B
258
       call pcondlsubr( S, T1V, ipar )
259
       residual = normfun( HUTI_NDIM, S, 1 ) / rhsnorm
260
    case (HUTI_PSEUDORESIDUAL)
261
       residual = normfun( HUTI_NDIM, R, 1 )
262
    case (HUTI_PRESID_SCALED_BYB)
263
       residual = normfun( HUTI_NDIM, R, 1 ) / rhsnorm
264
    case (HUTI_PRESID_SCALED_BYPRECB)
265
       residual = normfun( HUTI_NDIM, R, 1 ) / precrhsnorm
266
    case (HUTI_XDIFF_NORM)
267
       T1V = omega1 * R + omega2 * S + alpha * U
268
       residual = normfun( HUTI_NDIM, T1V, 1 )
269
    case (HUTI_USUPPLIED_STOPC)
270
       residual = stopcfun( X, B, R, ipar, dpar )
271
    case default
272
       call pcondrsubr( S, X, ipar )
273
       call matvecsubr( S, T1V, ipar )
274
       T1V = T1V - B
275
       call pcondlsubr( S, T1V, ipar )
276
       residual = normfun( HUTI_NDIM, S, 1 )
277
    end select
278

279
    !
280
    ! Print debugging output if desired
281
    !
282

283
    if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
284
       if ( mod(iter_count, HUTI_DBUGLVL) .eq. 0 ) then
285
          write (*, '(I8, E11.4)') iter_count, residual
286
       end if
287
    end if
288

289
    if ( residual .lt. HUTI_TOLERANCE ) then
290
       HUTI_INFO = HUTI_CONVERGENCE
291
       go to 1000
292
    end if
293

294
    IF( residual /= residual .OR. residual > HUTI_MAXTOLERANCE ) THEN
295
      HUTI_INFO = HUTI_DIVERGENCE
296
      GOTO 1000
297
    END IF
298
    
299
    U = U - omega1 * V - omega2 * W
300

301
    !
302
    ! Return back to the iteration loop (without initialization)
303
    !
304

305
    iter_count = iter_count + 1
306
    if ( iter_count .gt. HUTI_MAXIT ) then
307
       HUTI_INFO = HUTI_MAXITER
308
       go to 1000
309
    end if
310

311
    go to 300
312

313
    !
314
    ! This is where we exit last time (after enough iterations or breakdown)
315
    !
316

317
1000 continue
318
    if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
319
       write (*, '(I8, E11.4)') iter_count, residual
320
    end if
321

322
    HUTI_ITERS = iter_count
323
    return
324

325
    ! End of execution
326
    !*********************************************************************
327

328
  end subroutine  huti_sbicgstab_2solv
329

330
  !*************************************************************************
331

332

333
  !*************************************************************************
334
  !*************************************************************************
335
  ! Double precision version
336
  !*************************************************************************
337
  !*************************************************************************
338

339
  subroutine  huti_dbicgstab_2solv  ( ndim, wrkdim, xvec, rhsvec, &
340
       ipar, dpar, work, matvecsubr, pcondlsubr, &
341
       pcondrsubr, dotprodfun, normfun, stopcfun )
342

343
    use huti_interfaces
344
    implicit none
345

346
    procedure( mv_iface_d ), pointer :: matvecsubr 
347
    procedure( pc_iface_d ), pointer :: pcondlsubr 
348
    procedure( pc_iface_d ), pointer :: pcondrsubr 
349
    procedure( dotp_iface_d ), pointer :: dotprodfun 
350
    procedure( norm_iface_d ), pointer :: normfun 
351
    procedure( stopc_iface_d ), pointer :: stopcfun 
352

353
    ! Parameters
354

355
    integer :: ndim, wrkdim
356
    double precision, dimension(ndim) :: xvec, rhsvec
357
    integer, dimension(HUTI_IPAR_DFLTSIZE) :: ipar
358
    double precision, dimension(HUTI_DPAR_DFLTSIZE) :: dpar
359
    double precision, dimension(ndim,wrkdim) :: work
360

361
    ! Local variables
362

363
    double precision :: rho, oldrho, alpha, beta, omega1, omega2
364
    double precision :: tau, delta, myy
365
    integer :: iter_count
366

367
    double precision :: residual, rhsnorm, precrhsnorm
368

369
    !
370
    ! End of variable declarations
371
    !*********************************************************************
372

373
    !*********************************************************************
374
    ! The actual BiCGSTAB begins here (look the pseudo code in the
375
    ! "Templates..."-book, page 27)
376
    !
377
    ! First the initialization part
378
    !
379

380
    iter_count = 1
381

382
    ! The following applies for all matrix operations in this solver
383

384
    HUTI_EXTOP_MATTYPE = HUTI_MAT_NOTTRPSED
385

386
    ! Norms of right-hand side vector are used in convergence tests
387

388
    if ( HUTI_STOPC .eq. HUTI_TRESID_SCALED_BYB .or. & 
389
         HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYB ) then
390
       rhsnorm = normfun( HUTI_NDIM, B, 1 )
391
    end if
392
    if ( HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYPRECB ) then
393
       call pcondlsubr( T1V, B, ipar )
394
       precrhsnorm = normfun( HUTI_NDIM, T1V, 1 )
395
    end if
396

397
    ! Generate vector X if needed
398

399
    if ( HUTI_INITIALX .eq. HUTI_RANDOMX ) then
400
       call  huti_drandvec   ( X, ipar )
401
    else if ( HUTI_INITIALX .ne. HUTI_USERSUPPLIEDX ) then
402
       X = 1
403
    end if
404

405
    call pcondrsubr( U, X, ipar )
406
    call matvecsubr( U, R, ipar )
407
    U = B - R
408
    call pcondlsubr( R, U, ipar )
409
    RTLD = R
410
    U = 0
411
    oldrho = 1; omega2 = 1; alpha = 0
412

413
    !
414
    ! This is where the loop starts (that is we continue from here after
415
    ! the first iteration)
416
    !
417

418
300 continue
419

420
    oldrho = -omega2 * oldrho
421

422
    !
423
    ! This is the even BiCG step
424
    !
425

426
    rho = dotprodfun( HUTI_NDIM, RTLD, 1, R, 1 )
427
    if ( rho .eq. 0 ) then
428
       HUTI_INFO = HUTI_BICGSTAB_2_RHO
429
       go to 1000
430
    end if
431

432
    beta = ( rho * alpha ) / oldrho
433
    oldrho = rho
434
    U = R - beta * U
435

436
    call pcondrsubr( V, U, ipar )
437
    call matvecsubr( V, T1V, ipar )
438
    call pcondlsubr( V, T1V, ipar )
439

440
    alpha = oldrho / dotprodfun( HUTI_NDIM, RTLD, 1, V, 1 )
441
    R = R - alpha * V
442

443
    call pcondrsubr( S, R, ipar )
444
    call matvecsubr( S, T1V, ipar )
445
    call pcondlsubr( S, T1V, ipar )
446

447
    X = X + alpha * U
448

449
    !
450
    ! This is the odd BiCG step
451
    !
452

453
    rho = dotprodfun( HUTI_NDIM, RTLD, 1, S, 1 )
454
    if ( rho .eq. 0 ) then
455
       HUTI_INFO = HUTI_BICGSTAB_2_RHO
456
       go to 1000
457
    end if
458

459
    beta = ( rho * alpha ) / oldrho
460
    oldrho = rho
461
    V = S - beta * V
462

463
    call pcondrsubr( W, V, ipar )
464
    call matvecsubr( W, T1V, ipar )
465
    call pcondlsubr( W, T1V, ipar )
466

467
    alpha = oldrho / dotprodfun( HUTI_NDIM, RTLD, 1, W, 1 )
468
    U = R - beta * U
469
    R = R - alpha * V
470
    S = S - alpha * W
471

472
    call pcondrsubr( T, S, ipar )
473
    call matvecsubr( T, T1V, ipar )
474
    call pcondlsubr( T, T1V, ipar )
475

476
    !
477
    ! This is the GCR(2) part
478
    !
479

480
    omega1 = dotprodfun( HUTI_NDIM, R, 1, S, 1 )
481
    myy = dotprodfun( HUTI_NDIM, S, 1, S, 1 )
482
    delta = dotprodfun( HUTI_NDIM, S, 1, T, 1 )
483
    tau = dotprodfun( HUTI_NDIM, T, 1, T, 1 )
484
    omega2 = dotprodfun( HUTI_NDIM, R, 1, T, 1 )
485

486
    tau = tau - ( delta * delta ) / myy
487
    omega2 = ( omega2 - ( delta * omega1 ) / myy ) / tau
488
    omega1 = ( omega1 - delta * omega2 ) / myy
489

490
    X = X + omega1 * R + omega2 * S + alpha * U
491
    R = R - omega1 * S - omega2 * T
492

493
    !
494
    ! Check the convergence against selected stopping criterion
495
    !
496

497
    select case (HUTI_STOPC)
498
    case (HUTI_TRUERESIDUAL)
499
       call pcondrsubr( S, X, ipar )
500
       call matvecsubr( S, T1V, ipar )
501
       T1V = T1V - B
502
       call pcondlsubr( S, T1V, ipar )
503
       residual = normfun( HUTI_NDIM, S, 1 )
504
    case (HUTI_TRESID_SCALED_BYB)
505
       call pcondrsubr( S, X, ipar )
506
       call matvecsubr( S, T1V, ipar )
507
       T1V = T1V - B
508
       call pcondlsubr( S, T1V, ipar )
509
       residual = normfun( HUTI_NDIM, S, 1 ) / rhsnorm
510
    case (HUTI_PSEUDORESIDUAL)
511
       residual = normfun( HUTI_NDIM, R, 1 )
512
    case (HUTI_PRESID_SCALED_BYB)
513
       residual = normfun( HUTI_NDIM, R, 1 ) / rhsnorm
514
    case (HUTI_PRESID_SCALED_BYPRECB)
515
       residual = normfun( HUTI_NDIM, R, 1 ) / precrhsnorm
516
    case (HUTI_XDIFF_NORM)
517
       T1V = omega1 * R + omega2 * S + alpha * U
518
       residual = normfun( HUTI_NDIM, T1V, 1 )
519
    case (HUTI_USUPPLIED_STOPC)
520
       residual = stopcfun( X, B, R, ipar, dpar )
521
    case default
522
       call pcondrsubr( S, X, ipar )
523
       call matvecsubr( S, T1V, ipar )
524
       T1V = T1V - B
525
       call pcondlsubr( S, T1V, ipar )
526
       residual = normfun( HUTI_NDIM, S, 1 )
527
    end select
528

529
    !
530
    ! Print debugging output if desired
531
    !
532

533
    if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
534
       if ( mod(iter_count, HUTI_DBUGLVL) .eq. 0 ) then
535
          write (*, '(I8, E11.4)') iter_count, residual
536
       end if
537
    end if
538

539
    if ( residual .lt. HUTI_TOLERANCE ) then
540
       HUTI_INFO = HUTI_CONVERGENCE
541
       go to 1000
542
    end if
543

544
    IF( residual /= residual .OR. residual > HUTI_MAXTOLERANCE ) THEN
545
      HUTI_INFO = HUTI_DIVERGENCE
546
      GOTO 1000
547
    END IF
548
    
549
    U = U - omega1 * V - omega2 * W
550

551
    !
552
    ! Return back to the iteration loop (without initialization)
553
    !
554

555
    iter_count = iter_count + 1
556
    if ( iter_count .gt. HUTI_MAXIT ) then
557
       HUTI_INFO = HUTI_MAXITER
558
       go to 1000
559
    end if
560

561
    go to 300
562

563
    !
564
    ! This is where we exit last time (after enough iterations or breakdown)
565
    !
566

567
1000 continue
568
    if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
569
       write (*, '(I8, E11.4)') iter_count, residual
570
    end if
571

572
    HUTI_ITERS = iter_count
573
    return
574

575
    ! End of execution
576
    !*********************************************************************
577

578
  end subroutine  huti_dbicgstab_2solv
579

580
  !*************************************************************************
581

582

583
  !*************************************************************************
584
  !*************************************************************************
585
  ! Complex version
586
  !*************************************************************************
587
  !*************************************************************************
588

589
  subroutine  huti_cbicgstab_2solv  ( ndim, wrkdim, xvec, rhsvec, &
590
       ipar, dpar, work, matvecsubr, pcondlsubr, &
591
       pcondrsubr, dotprodfun, normfun, stopcfun )
592

593
    use huti_interfaces
594
    implicit none
595

596
    procedure( mv_iface_c ), pointer :: matvecsubr 
597
    procedure( pc_iface_c ), pointer :: pcondlsubr 
598
    procedure( pc_iface_c ), pointer :: pcondrsubr 
599
    procedure( dotp_iface_c ), pointer :: dotprodfun 
600
    procedure( norm_iface_c ), pointer :: normfun 
601
    procedure( stopc_iface_c ), pointer :: stopcfun 
602

603
    ! Parameters
604

605
    integer :: ndim, wrkdim
606
    complex, dimension(ndim) :: xvec, rhsvec
607
    integer, dimension(HUTI_IPAR_DFLTSIZE) :: ipar
608
    double precision, dimension(HUTI_DPAR_DFLTSIZE) :: dpar
609
    complex, dimension(ndim,wrkdim) :: work
610

611
    ! Local variables
612

613
    complex :: rho, oldrho, alpha, beta, omega1, omega2
614
    complex :: tau, delta, myy
615
    integer :: iter_count
616

617
    real :: residual, rhsnorm, precrhsnorm
618

619
    !
620
    ! End of variable declarations
621
    !*********************************************************************
622

623
    !*********************************************************************
624
    ! The actual BiCGSTAB begins here (look the pseudo code in the
625
    ! "Templates..."-book, page 27)
626
    !
627
    ! First the initialization part
628
    !
629

630
    iter_count = 1
631

632
    ! The following applies for all matrix operations in this solver
633

634
    HUTI_EXTOP_MATTYPE = HUTI_MAT_NOTTRPSED
635

636
    ! Norms of right-hand side vector are used in convergence tests
637

638
    if ( HUTI_STOPC .eq. HUTI_TRESID_SCALED_BYB .or. & 
639
         HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYB ) then
640
       rhsnorm = normfun( HUTI_NDIM, B, 1 )
641
    end if
642
    if ( HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYPRECB ) then
643
       call pcondlsubr( T1V, B, ipar )
644
       precrhsnorm = normfun( HUTI_NDIM, T1V, 1 )
645
    end if
646

647
    ! Generate vector X if needed
648

649
    if ( HUTI_INITIALX .eq. HUTI_RANDOMX ) then
650
       call  huti_crandvec   ( X, ipar )
651
    else if ( HUTI_INITIALX .ne. HUTI_USERSUPPLIEDX ) then
652
       X = 1
653
    end if
654

655
    call pcondrsubr( U, X, ipar )
656
    call matvecsubr( U, R, ipar )
657
    U = B - R
658
    call pcondlsubr( R, U, ipar )
659
    RTLD = R
660
    U = 0
661
    oldrho = 1; omega2 = 1; alpha = 0
662

663
    !
664
    ! This is where the loop starts (that is we continue from here after
665
    ! the first iteration)
666
    !
667

668
300 continue
669

670
    oldrho = -omega2 * oldrho
671

672
    !
673
    ! This is the even BiCG step
674
    !
675

676
    rho = dotprodfun( HUTI_NDIM, RTLD, 1, R, 1 )
677
    if ( rho .eq. 0 ) then
678
       HUTI_INFO = HUTI_BICGSTAB_2_RHO
679
       go to 1000
680
    end if
681

682
    beta = ( rho * alpha ) / oldrho
683
    oldrho = rho
684
    U = R - beta * U
685

686
    call pcondrsubr( V, U, ipar )
687
    call matvecsubr( V, T1V, ipar )
688
    call pcondlsubr( V, T1V, ipar )
689

690
    alpha = oldrho / dotprodfun( HUTI_NDIM, RTLD, 1, V, 1 )
691
    R = R - alpha * V
692

693
    call pcondrsubr( S, R, ipar )
694
    call matvecsubr( S, T1V, ipar )
695
    call pcondlsubr( S, T1V, ipar )
696

697
    X = X + alpha * U
698

699
    !
700
    ! This is the odd BiCG step
701
    !
702

703
    rho = dotprodfun( HUTI_NDIM, RTLD, 1, S, 1 )
704
    if ( rho .eq. 0 ) then
705
       HUTI_INFO = HUTI_BICGSTAB_2_RHO
706
       go to 1000
707
    end if
708

709
    beta = ( rho * alpha ) / oldrho
710
    oldrho = rho
711
    V = S - beta * V
712

713
    call pcondrsubr( W, V, ipar )
714
    call matvecsubr( W, T1V, ipar )
715
    call pcondlsubr( W, T1V, ipar )
716

717
    alpha = oldrho / dotprodfun( HUTI_NDIM, RTLD, 1, W, 1 )
718
    U = R - beta * U
719
    R = R - alpha * V
720
    S = S - alpha * W
721

722
    call pcondrsubr( T, S, ipar )
723
    call matvecsubr( T, T1V, ipar )
724
    call pcondlsubr( T, T1V, ipar )
725

726
    !
727
    ! This is the GCR(2) part
728
    !
729

730
    omega1 = dotprodfun( HUTI_NDIM, R, 1, S, 1 )
731
    myy = dotprodfun( HUTI_NDIM, S, 1, S, 1 )
732
    delta = dotprodfun( HUTI_NDIM, S, 1, T, 1 )
733
    tau = dotprodfun( HUTI_NDIM, T, 1, T, 1 )
734
    omega2 = dotprodfun( HUTI_NDIM, R, 1, T, 1 )
735

736
    tau = tau - ( delta * delta ) / myy
737
    omega2 = ( omega2 - ( delta * omega1 ) / myy ) / tau
738
    omega1 = ( omega1 - delta * omega2 ) / myy
739

740
    X = X + omega1 * R + omega2 * S + alpha * U
741
    R = R - omega1 * S - omega2 * T
742

743
    !
744
    ! Check the convergence against selected stopping criterion
745
    !
746

747
    select case (HUTI_STOPC)
748
    case (HUTI_TRUERESIDUAL)
749
       call pcondrsubr( S, X, ipar )
750
       call matvecsubr( S, T1V, ipar )
751
       T1V = T1V - B
752
       call pcondlsubr( S, T1V, ipar )
753
       residual = normfun( HUTI_NDIM, S, 1 )
754
    case (HUTI_TRESID_SCALED_BYB)
755
       call pcondrsubr( S, X, ipar )
756
       call matvecsubr( S, T1V, ipar )
757
       T1V = T1V - B
758
       call pcondlsubr( S, T1V, ipar )
759
       residual = normfun( HUTI_NDIM, S, 1 ) / rhsnorm
760
    case (HUTI_PSEUDORESIDUAL)
761
       residual = normfun( HUTI_NDIM, R, 1 )
762
    case (HUTI_PRESID_SCALED_BYB)
763
       residual = normfun( HUTI_NDIM, R, 1 ) / rhsnorm
764
    case (HUTI_PRESID_SCALED_BYPRECB)
765
       residual = normfun( HUTI_NDIM, R, 1 ) / precrhsnorm
766
    case (HUTI_XDIFF_NORM)
767
       T1V = omega1 * R + omega2 * S + alpha * U
768
       residual = normfun( HUTI_NDIM, T1V, 1 )
769
    case (HUTI_USUPPLIED_STOPC)
770
       residual = stopcfun( X, B, R, ipar, dpar )
771
    case default
772
       call pcondrsubr( S, X, ipar )
773
       call matvecsubr( S, T1V, ipar )
774
       T1V = T1V - B
775
       call pcondlsubr( S, T1V, ipar )
776
       residual = normfun( HUTI_NDIM, S, 1 )
777
    end select
778

779
    !
780
    ! Print debugging output if desired
781
    !
782

783
    if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
784
       if ( mod(iter_count, HUTI_DBUGLVL) .eq. 0 ) then
785
          write (*, '(I8, E11.4)') iter_count, residual
786
       end if
787
    end if
788

789
    if ( residual .lt. HUTI_TOLERANCE ) then
790
       HUTI_INFO = HUTI_CONVERGENCE
791
       go to 1000
792
    end if
793

794
    IF( residual /= residual .OR. residual > HUTI_MAXTOLERANCE ) THEN
795
      HUTI_INFO = HUTI_DIVERGENCE
796
      GOTO 1000
797
    END IF
798

799

800
    
801
    U = U - omega1 * V - omega2 * W
802

803
    !
804
    ! Return back to the iteration loop (without initialization)
805
    !
806

807
    iter_count = iter_count + 1
808
    if ( iter_count .gt. HUTI_MAXIT ) then
809
       HUTI_INFO = HUTI_MAXITER
810
       go to 1000
811
    end if
812

813
    go to 300
814

815
    !
816
    ! This is where we exit last time (after enough iterations or breakdown)
817
    !
818

819
1000 continue
820
    if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
821
       write (*, '(I8, E11.4)') iter_count, residual
822
    end if
823

824
    HUTI_ITERS = iter_count
825
    return
826

827
    ! End of execution
828
    !*********************************************************************
829

830
  end subroutine  huti_cbicgstab_2solv
831

832
  !*************************************************************************
833

834

835
  !*************************************************************************
836
  !*************************************************************************
837
  ! Double complex version
838
  !*************************************************************************
839
  !*************************************************************************
840

841
  subroutine  huti_zbicgstab_2solv  ( ndim, wrkdim, xvec, rhsvec, &
842
       ipar, dpar, work, matvecsubr, pcondlsubr, &
843
       pcondrsubr, dotprodfun, normfun, stopcfun )
844

845
    use huti_interfaces
846
    implicit none
847

848
    procedure( mv_iface_z ), pointer :: matvecsubr 
849
    procedure( pc_iface_z ), pointer :: pcondlsubr
850
    procedure( pc_iface_z ), pointer :: pcondrsubr 
851
    procedure( dotp_iface_z ), pointer :: dotprodfun 
852
    procedure( norm_iface_z ), pointer :: normfun 
853
    procedure( stopc_iface_z ), pointer :: stopcfun 
854

855
    ! Parameters
856

857
    integer :: ndim, wrkdim
858
    double complex, dimension(ndim) :: xvec, rhsvec
859
    integer, dimension(HUTI_IPAR_DFLTSIZE) :: ipar
860
    double precision, dimension(HUTI_DPAR_DFLTSIZE) :: dpar
861
    double complex, dimension(ndim,wrkdim) :: work
862

863
    ! Local variables
864

865
    double complex :: rho, oldrho, alpha, beta, omega1, omega2
866
    double complex :: tau, delta, myy
867
    integer :: iter_count
868

869
    double precision :: residual, rhsnorm, precrhsnorm
870

871
    !
872
    ! End of variable declarations
873
    !*********************************************************************
874

875
    !*********************************************************************
876
    ! The actual BiCGSTAB begins here (look the pseudo code in the
877
    ! "Templates..."-book, page 27)
878
    !
879
    ! First the initialization part
880
    !
881

882
    iter_count = 1
883

884
    ! The following applies for all matrix operations in this solver
885

886
    HUTI_EXTOP_MATTYPE = HUTI_MAT_NOTTRPSED
887

888
    ! Norms of right-hand side vector are used in convergence tests
889

890
    if ( HUTI_STOPC .eq. HUTI_TRESID_SCALED_BYB .or. & 
891
         HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYB ) then
892
       rhsnorm = normfun( HUTI_NDIM, B, 1 )
893
    end if
894
    if ( HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYPRECB ) then
895
       call pcondlsubr( T1V, B, ipar )
896
       precrhsnorm = normfun( HUTI_NDIM, T1V, 1 )
897
    end if
898

899
    ! Generate vector X if needed
900

901
    if ( HUTI_INITIALX .eq. HUTI_RANDOMX ) then
902
       call  huti_zrandvec   ( X, ipar )
903
    else if ( HUTI_INITIALX .ne. HUTI_USERSUPPLIEDX ) then
904
       X = 1
905
    end if
906

907
    call pcondrsubr( U, X, ipar )
908
    call matvecsubr( U, R, ipar )
909
    U = B - R
910
    call pcondlsubr( R, U, ipar )
911
    RTLD = R
912
    U = 0
913
    oldrho = 1; omega2 = 1; alpha = 0
914

915
    !
916
    ! This is where the loop starts (that is we continue from here after
917
    ! the first iteration)
918
    !
919

920
300 continue
921

922
    oldrho = -omega2 * oldrho
923

924
    !
925
    ! This is the even BiCG step
926
    !
927

928
    rho = dotprodfun( HUTI_NDIM, RTLD, 1, R, 1 )
929
    if ( rho .eq. 0 ) then
930
       HUTI_INFO = HUTI_BICGSTAB_2_RHO
931
       go to 1000
932
    end if
933

934
    beta = ( rho * alpha ) / oldrho
935
    oldrho = rho
936
    U = R - beta * U
937

938
    call pcondrsubr( V, U, ipar )
939
    call matvecsubr( V, T1V, ipar )
940
    call pcondlsubr( V, T1V, ipar )
941

942
    alpha = oldrho / dotprodfun( HUTI_NDIM, RTLD, 1, V, 1 )
943
    R = R - alpha * V
944

945
    call pcondrsubr( S, R, ipar )
946
    call matvecsubr( S, T1V, ipar )
947
    call pcondlsubr( S, T1V, ipar )
948

949
    X = X + alpha * U
950

951
    !
952
    ! This is the odd BiCG step
953
    !
954

955
    rho = dotprodfun( HUTI_NDIM, RTLD, 1, S, 1 )
956
    if ( rho .eq. 0 ) then
957
       HUTI_INFO = HUTI_BICGSTAB_2_RHO
958
       go to 1000
959
    end if
960

961
    beta = ( rho * alpha ) / oldrho
962
    oldrho = rho
963
    V = S - beta * V
964

965
    call pcondrsubr( W, V, ipar )
966
    call matvecsubr( W, T1V, ipar )
967
    call pcondlsubr( W, T1V, ipar )
968

969
    alpha = oldrho / dotprodfun( HUTI_NDIM, RTLD, 1, W, 1 )
970
    U = R - beta * U
971
    R = R - alpha * V
972
    S = S - alpha * W
973

974
    call pcondrsubr( T, S, ipar )
975
    call matvecsubr( T, T1V, ipar )
976
    call pcondlsubr( T, T1V, ipar )
977

978
    !
979
    ! This is the GCR(2) part
980
    !
981

982
    omega1 = dotprodfun( HUTI_NDIM, R, 1, S, 1 )
983
    myy = dotprodfun( HUTI_NDIM, S, 1, S, 1 )
984
    delta = dotprodfun( HUTI_NDIM, S, 1, T, 1 )
985
    tau = dotprodfun( HUTI_NDIM, T, 1, T, 1 )
986
    omega2 = dotprodfun( HUTI_NDIM, R, 1, T, 1 )
987

988
    tau = tau - ( delta * delta ) / myy
989
    omega2 = ( omega2 - ( delta * omega1 ) / myy ) / tau
990
    omega1 = ( omega1 - delta * omega2 ) / myy
991

992
    X = X + omega1 * R + omega2 * S + alpha * U
993
    R = R - omega1 * S - omega2 * T
994

995
    !
996
    ! Check the convergence against selected stopping criterion
997
    !
998

999
    select case (HUTI_STOPC)
1000
    case (HUTI_TRUERESIDUAL)
1001
       call pcondrsubr( S, X, ipar )
1002
       call matvecsubr( S, T1V, ipar )
1003
       T1V = T1V - B
1004
       call pcondlsubr( S, T1V, ipar )
1005
       residual = normfun( HUTI_NDIM, S, 1 )
1006
    case (HUTI_TRESID_SCALED_BYB)
1007
       call pcondrsubr( S, X, ipar )
1008
       call matvecsubr( S, T1V, ipar )
1009
       T1V = T1V - B
1010
       call pcondlsubr( S, T1V, ipar )
1011
       residual = normfun( HUTI_NDIM, S, 1 ) / rhsnorm
1012
    case (HUTI_PSEUDORESIDUAL)
1013
       residual = normfun( HUTI_NDIM, R, 1 )
1014
    case (HUTI_PRESID_SCALED_BYB)
1015
       residual = normfun( HUTI_NDIM, R, 1 ) / rhsnorm
1016
    case (HUTI_PRESID_SCALED_BYPRECB)
1017
       residual = normfun( HUTI_NDIM, R, 1 ) / precrhsnorm
1018
    case (HUTI_XDIFF_NORM)
1019
       T1V = omega1 * R + omega2 * S + alpha * U
1020
       residual = normfun( HUTI_NDIM, T1V, 1 )
1021
    case (HUTI_USUPPLIED_STOPC)
1022
       residual = stopcfun( X, B, R, ipar, dpar )
1023
    case default
1024
       call pcondrsubr( S, X, ipar )
1025
       call matvecsubr( S, T1V, ipar )
1026
       T1V = T1V - B
1027
       call pcondlsubr( S, T1V, ipar )
1028
       residual = normfun( HUTI_NDIM, S, 1 )
1029
    end select
1030

1031
    !
1032
    ! Print debugging output if desired
1033
    !
1034

1035
    if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
1036
       if ( mod(iter_count, HUTI_DBUGLVL) .eq. 0 ) then
1037
          write (*, '(I8, E11.4)') iter_count, residual
1038
       end if
1039
    end if
1040

1041
    if ( residual .lt. HUTI_TOLERANCE ) then
1042
       HUTI_INFO = HUTI_CONVERGENCE
1043
       go to 1000
1044
    end if
1045

1046
    IF( residual /= residual .OR. residual > HUTI_MAXTOLERANCE ) THEN
1047
      HUTI_INFO = HUTI_DIVERGENCE
1048
      GOTO 1000
1049
    END IF
1050
    
1051
    U = U - omega1 * V - omega2 * W
1052

1053
    !
1054
    ! Return back to the iteration loop (without initialization)
1055
    !
1056

1057
    iter_count = iter_count + 1
1058
    if ( iter_count .gt. HUTI_MAXIT ) then
1059
       HUTI_INFO = HUTI_MAXITER
1060
       go to 1000
1061
    end if
1062

1063
    go to 300
1064

1065
    !
1066
    ! This is where we exit last time (after enough iterations or breakdown)
1067
    !
1068

1069
1000 continue
1070
    if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
1071
       write (*, '(I8, E11.4)') iter_count, residual
1072
    end if
1073

1074
    HUTI_ITERS = iter_count
1075
    return
1076

1077
    ! End of execution
1078
    !*********************************************************************
1079

1080
  end subroutine  huti_zbicgstab_2solv
1081

1082
  !*************************************************************************
1083

1084
end module huti_bicgstab_2
1085

1086
Product

Resources

Company
