CoCalc -- huti_bicgstab.F90

GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/fhutiter/src/huti_bicgstab.F90
³²⁰⁶ views
1
module huti_bicgstab
2
  use huti_aux
3
  implicit none
4

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

28
  !
29
  ! Subroutine to implement BiConjugate Gradient Stabilised iteration
30
  !
31
  ! $Id: huti_bicgstab.src,v 1.1.1.1 2005/04/15 10:31:18 vierinen Exp $
32

33

34
#include  "huti_fdefs.h" 
35

36
  !
37
  !*************************************************************************
38
  ! Definitions to make the code more understandable and to make it look
39
  ! like the pseudo code
40
  !
41

42
#define  X  xvec 
43
#define  B  rhsvec 
44

45
#define  RTLD  work(:,1) 
46
#define  RTLD_ind  1 
47
#define  P  work(:,2) 
48
#define  P_ind  2 
49
#define  T1V  work(:,3) 
50
#define  T1V_ind  3 
51
#define  V  work(:,4) 
52
#define  V_ind  4 
53
#define  S  work(:,5) 
54
#define  S_ind  5 
55
#define  T2V  work(:,6) 
56
#define  T2V_ind  6 
57
#define  T  work(:,7) 
58
#define  T_ind  7 
59
#define  R  work(:,8) 
60
#define  R_ind  8 
61

62

63
contains
64

65

66
  !*************************************************************************
67
  !*************************************************************************
68
  ! Single precision version
69
  !*************************************************************************
70
  !*************************************************************************
71

72
  subroutine  huti_sbicgstabsolv  ( ndim, wrkdim, xvec, rhsvec, &
73
       ipar, dpar, work, matvecsubr, pcondlsubr, &
74
       pcondrsubr, dotprodfun, normfun, stopcfun )
75

76
    use huti_interfaces
77
    implicit none
78

79
    procedure( mv_iface_s ), pointer :: matvecsubr
80
    procedure( pc_iface_s ), pointer :: pcondlsubr
81
    procedure( pc_iface_s ), pointer :: pcondrsubr
82
    procedure( dotp_iface_s ), pointer :: dotprodfun
83
    procedure( norm_iface_s ), pointer :: normfun
84
    procedure( stopc_iface_s ), pointer :: stopcfun
85

86
    ! Parameters
87

88
    integer :: ndim, wrkdim
89
    real, dimension(ndim) :: xvec, rhsvec
90
    integer, dimension(HUTI_IPAR_DFLTSIZE) :: ipar
91
    double precision, dimension(HUTI_DPAR_DFLTSIZE) :: dpar
92
    real, dimension(ndim,wrkdim) :: work
93

94
    ! Local variables
95

96
    real :: rho, oldrho, alpha, beta, omega
97
    integer :: iter_count
98

99
    real :: residual, rhsnorm, precrhsnorm
100

101
    !
102
    ! End of variable declarations
103
    !*********************************************************************
104

105
    !*********************************************************************
106
    ! The actual BiCGSTAB begins here (look the pseudo code in the
107
    ! "Templates..."-book, page 27)
108
    !
109
    ! First the initialization part
110
    !
111

112
    iter_count = 1
113

114
    ! The following applies for all matrix operations in this solver
115

116
    HUTI_EXTOP_MATTYPE = HUTI_MAT_NOTTRPSED
117

118
    ! Norms of right-hand side vector are used in convergence tests
119

120
    if ( HUTI_STOPC .eq. HUTI_TRESID_SCALED_BYB .or. & 
121
         HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYB ) then
122
       rhsnorm = normfun( HUTI_NDIM, B, 1 )
123
    end if
124
    if ( HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYPRECB ) then
125
       call pcondlsubr( P, B, ipar )
126
       precrhsnorm = normfun( HUTI_NDIM, P, 1 )
127
    end if
128

129
    ! Generate vector X if needed
130

131
    if ( HUTI_INITIALX .eq. HUTI_RANDOMX ) then
132
       call  huti_srandvec   ( X, ipar )
133
    else if ( HUTI_INITIALX .ne. HUTI_USERSUPPLIEDX ) then
134
       X = 1
135
    end if
136

137
    call matvecsubr( X, R, ipar )
138
    R = B - R
139
    RTLD = R
140
    P = 0; V = 0
141
    oldrho = 1; omega = 1; alpha = 0
142

143
    !
144
    ! This is where the loop starts (that is we continue from here after
145
    ! the first iteration)
146
    !
147

148
300 continue
149

150
    rho = dotprodfun( HUTI_NDIM, RTLD, 1, R, 1 )
151
    if ( rho .eq. 0 ) then
152
       HUTI_INFO = HUTI_BICGSTAB_RHO
153
       go to 1000
154
    end if
155

156
    beta = ( rho * alpha ) / ( oldrho * omega )
157
    P = R + beta * ( P - omega * V )
158

159
    call pcondlsubr( V, P, ipar )
160
    call pcondrsubr( T1V, V, ipar )
161
    call matvecsubr( T1V, V, ipar )
162

163
    alpha = rho / dotprodfun( HUTI_NDIM, RTLD, 1, V, 1 )
164
    S = R - alpha * V
165

166
    residual = normfun( HUTI_NDIM, S, 1 )
167
    if ( residual .lt. HUTI_EPSILON ) then
168
       X = X + alpha * T1V
169
       !HUTI_INFO = HUTI_BICGSTAB_SNORM
170
       HUTI_INFO = HUTI_CONVERGENCE
171
       go to 1000
172
    end if
173

174
    call pcondlsubr( T, S, ipar )
175
    call pcondrsubr( T2V, T, ipar )
176
    call matvecsubr( T2V, T, ipar )
177

178
    omega = ( dotprodfun( HUTI_NDIM, T, 1, S, 1 ) ) / &
179
         ( dotprodfun( HUTI_NDIM, T, 1, T, 1 ) )
180
    X = X + alpha * T1V + omega * T2V
181
    R = S - omega * T
182

183
    !
184
    ! Check the convergence against selected stopping criterion
185
    !
186

187
    select case (HUTI_STOPC)
188
    case (HUTI_TRUERESIDUAL)
189
       call matvecsubr( X, T2V, ipar )
190
       T1V = T2V - B
191
       residual = normfun( HUTI_NDIM, T1V, 1 )
192
    case (HUTI_TRESID_SCALED_BYB)
193
       call matvecsubr( X, T2V, ipar )
194
       T1V = T2V - B
195
       residual = normfun( HUTI_NDIM, T1V, 1 ) / rhsnorm
196
    case (HUTI_PSEUDORESIDUAL)
197
       residual = normfun( HUTI_NDIM, R, 1 )
198
    case (HUTI_PRESID_SCALED_BYB)
199
       residual = normfun( HUTI_NDIM, R, 1 ) / rhsnorm
200
    case (HUTI_PRESID_SCALED_BYPRECB)
201
       residual = normfun( HUTI_NDIM, R, 1 ) / precrhsnorm
202
    case (HUTI_XDIFF_NORM)
203
       T1V = alpha * T1V + omega * T2V
204
       residual = normfun( HUTI_NDIM, T1V, 1 )
205
    case (HUTI_USUPPLIED_STOPC)
206
       residual = stopcfun( X, B, R, ipar, dpar )
207
    case default
208
       call matvecsubr( X, T2V, ipar )
209
       T1V = T2V - B
210
       residual = normfun( HUTI_NDIM, T1V, 1 )
211
    end select
212

213
    !
214
    ! Print debugging output if desired
215
    !
216

217
    if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
218
       if ( mod(iter_count, HUTI_DBUGLVL) .eq. 0 ) then
219
          write (*, '(I8, E11.4)') iter_count, residual
220
       end if
221
    end if
222

223
    if ( residual .lt. HUTI_TOLERANCE ) then
224
       HUTI_INFO = HUTI_CONVERGENCE
225
       go to 1000
226
    end if
227
    
228
    IF( residual /= residual .OR. residual > HUTI_MAXTOLERANCE ) THEN
229
      HUTI_INFO = HUTI_DIVERGENCE
230
      GOTO 1000
231
    END IF
232
    
233
    if ( omega .eq. 0 ) then
234
       HUTI_INFO = HUTI_BICGSTAB_OMEGA
235
       go to 1000
236
    end if
237

238
    oldrho = rho
239

240
    !
241
    ! Return back to the iteration loop (without initialization)
242
    !
243

244
    iter_count = iter_count + 1
245
    if ( iter_count .gt. HUTI_MAXIT ) then
246
       HUTI_INFO = HUTI_MAXITER
247
       go to 1000
248
    end if
249

250
    go to 300
251

252
    !
253
    ! This is where we exit last time (after enough iterations or breakdown)
254
    !
255

256
1000 continue
257
    if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
258
       write (*, '(I8, E11.4)') iter_count, residual
259
    end if
260

261
    HUTI_ITERS = iter_count
262
    return
263

264
    ! End of execution
265
    !*********************************************************************
266

267
  end subroutine  huti_sbicgstabsolv
268

269
  !*************************************************************************
270

271

272

273
  !*************************************************************************
274
  !*************************************************************************
275
  ! Double precision version
276
  !*************************************************************************
277
  !*************************************************************************
278

279
  subroutine  huti_dbicgstabsolv  ( ndim, wrkdim, xvec, rhsvec, &
280
       ipar, dpar, work, matvecsubr, pcondlsubr, &
281
       pcondrsubr, dotprodfun, normfun, stopcfun )
282

283
    use huti_interfaces
284
    implicit none
285

286
    procedure( mv_iface_d ), pointer :: matvecsubr
287
    procedure( pc_iface_d ), pointer :: pcondlsubr
288
    procedure( pc_iface_d ), pointer :: pcondrsubr
289
    procedure( dotp_iface_d ), pointer :: dotprodfun
290
    procedure( norm_iface_d ), pointer :: normfun
291
    procedure( stopc_iface_d ), pointer :: stopcfun
292

293
    ! Parameters
294

295
    integer :: ndim, wrkdim
296
    double precision, dimension(ndim) :: xvec, rhsvec
297
    integer, dimension(HUTI_IPAR_DFLTSIZE) :: ipar
298
    double precision, dimension(HUTI_DPAR_DFLTSIZE) :: dpar
299
    double precision, dimension(ndim,wrkdim) :: work
300

301
    ! Local variables
302

303
    double precision :: rho, oldrho, alpha, beta, omega
304
    integer :: iter_count, i
305

306
    double precision :: residual, rhsnorm, precrhsnorm
307

308
    !
309
    ! End of variable declarations
310
    !*********************************************************************
311

312
    !*********************************************************************
313
    ! The actual BiCGSTAB begins here (look the pseudo code in the
314
    ! "Templates..."-book, page 27)
315
    !
316
    ! First the initialization part
317
    !
318

319
    iter_count = 1
320

321
    ! The following applies for all matrix operations in this solver
322

323
    HUTI_EXTOP_MATTYPE = HUTI_MAT_NOTTRPSED
324

325
    ! Norms of right-hand side vector are used in convergence tests
326

327
    if ( HUTI_STOPC .eq. HUTI_TRESID_SCALED_BYB .or. & 
328
         HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYB ) then
329
       rhsnorm = normfun( HUTI_NDIM, B, 1 )
330
    end if
331
    if ( HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYPRECB ) then
332
       call pcondlsubr( P, B, ipar )
333
       precrhsnorm = normfun( HUTI_NDIM, P, 1 )
334
    end if
335

336
    ! Generate vector X if needed
337

338
    if ( HUTI_INITIALX .eq. HUTI_RANDOMX ) then
339
       call  huti_drandvec   ( X, ipar )
340
    else if ( HUTI_INITIALX .ne. HUTI_USERSUPPLIEDX ) then
341
#ifdef _OPENMP
342
       !$OMP PARALLEL DO
343
       do i=1,ndim
344
          xvec(i) = 1
345
       end do
346
       !$OMP END PARALLEL DO
347
#else
348
       X = 1
349
#endif
350
    end if
351

352
    call matvecsubr( X, R, ipar )
353
#ifdef _OPENMP
354
    !$OMP PARALLEL PRIVATE(i)
355
    !$OMP DO
356
    do i=1,ndim
357
       work(i,R_ind) = rhsvec(i) - work(i,R_ind)
358
       work(i,RTLD_ind) = work(i,R_ind)
359
    end do
360
    !$OMP END DO NOWAIT
361
    !$OMP DO
362
    do i=1,ndim
363
       work(i,P_ind) = 0
364
       work(i,V_ind) = 0
365
    end do
366
    !$OMP END DO NOWAIT
367
    !$OMP END PARALLEL
368
#else
369
    R = B - R
370
    RTLD = R
371
    P = 0; V = 0
372
#endif
373
    oldrho = 1; omega = 1; alpha = 0
374

375
    !
376
    ! This is where the loop starts (that is we continue from here after
377
    ! the first iteration)
378
    !
379

380
300 continue
381

382
    rho = dotprodfun( HUTI_NDIM, RTLD, 1, R, 1 )
383
    if ( rho .eq. 0 ) then
384
       HUTI_INFO = HUTI_BICGSTAB_RHO
385
       go to 1000
386
    end if
387

388
    beta = ( rho * alpha ) / ( oldrho * omega )
389

390
#ifdef _OPENMP
391
    !$OMP PARALLEL DO
392
    do i=1,ndim
393
       work(i,P_ind) = work(i,R_ind) + beta * ( work(i,P_ind) - omega * work(i,V_ind) )
394
    end do
395
    !$OMP END PARALLEL DO
396
#else
397
    P = R + beta * ( P - omega * V )
398
#endif
399

400
    call pcondlsubr( V, P, ipar )
401
    call pcondrsubr( T1V, V, ipar )
402
    call matvecsubr( T1V, V, ipar )
403

404
    alpha = rho / dotprodfun( HUTI_NDIM, RTLD, 1, V, 1 )
405
#ifdef _OPENMP
406
    !$OMP PARALLEL DO
407
    do i=1,ndim
408
       work(i,S_ind) = work(i,R_ind) - alpha * work(i,V_ind)
409
    end do
410
    !$OMP END PARALLEL DO
411
#else
412
    S = R - alpha * V
413
#endif
414

415
    residual = normfun( HUTI_NDIM, S, 1 )
416
    if ( residual .lt. HUTI_EPSILON ) then
417
#ifdef _OPENMP
418
    !$OMP PARALLEL DO
419
    do i=1,ndim
420
       xvec(i) = xvec(i) + alpha * work(i,T1V_ind)
421
    end do
422
    !$OMP END PARALLEL DO
423
#else
424
       X = X + alpha * T1V
425
#endif
426
       HUTI_INFO = HUTI_CONVERGENCE
427
       !HUTI_INFO = HUTI_BICGSTAB_SNORM
428
       go to 1000
429
    end if
430

431
    call pcondlsubr( T, S, ipar )
432
    call pcondrsubr( T2V, T, ipar )
433
    call matvecsubr( T2V, T, ipar )
434

435
    omega = ( dotprodfun( HUTI_NDIM, T, 1, S, 1 ) ) / &
436
         ( dotprodfun( HUTI_NDIM, T, 1, T, 1 ) )
437

438
#ifdef _OPENMP
439
    !$OMP PARALLEL DO
440
    do i=1,ndim
441
       xvec(i) = xvec(i) + alpha * work(i,T1V_ind) + omega * work(i,T2V_ind)
442
       work(i,R_ind) = work(i,S_ind) - omega * work(i,T_ind)
443
    end do
444
    !$OMP END PARALLEL DO
445
#else
446
    X = X + alpha * T1V + omega * T2V
447
    R = S - omega * T
448
#endif
449

450
    !
451
    ! Check the convergence against selected stopping criterion
452
    !
453

454
    select case (HUTI_STOPC)
455
    case (HUTI_TRUERESIDUAL)
456
       call matvecsubr( X, T2V, ipar )
457
#ifdef _OPENMP
458
       !$OMP PARALLEL DO
459
       do i=1,ndim
460
          work(i,T1V_ind) = work(i,T2V_ind) - rhsvec(i)
461
       end do
462
       !$OMP END PARALLEL DO
463
#else
464
       T1V = T2V - B
465
#endif
466
       residual = normfun( HUTI_NDIM, T1V, 1 )
467
    case (HUTI_TRESID_SCALED_BYB)
468
       call matvecsubr( X, T2V, ipar )
469
#ifdef _OPENMP
470
       !$OMP PARALLEL DO
471
       do i=1,ndim
472
          work(i,T1V_ind) = work(i,T2V_ind) - rhsvec(i)
473
       end do
474
       !$OMP END PARALLEL DO
475
#else
476
       T1V = T2V - B
477
#endif
478
       residual = normfun( HUTI_NDIM, T1V, 1 ) / rhsnorm
479
    case (HUTI_PSEUDORESIDUAL)
480
       residual = normfun( HUTI_NDIM, R, 1 )
481
    case (HUTI_PRESID_SCALED_BYB)
482
       residual = normfun( HUTI_NDIM, R, 1 ) / rhsnorm
483
    case (HUTI_PRESID_SCALED_BYPRECB)
484
       residual = normfun( HUTI_NDIM, R, 1 ) / precrhsnorm
485
    case (HUTI_XDIFF_NORM)
486
#ifdef _OPENMP
487
       !$OMP PARALLEL DO
488
       do i=1,ndim
489
          work(i,T1V_ind) = alpha * work(i,T1V_ind) + omega * work(i,T2V_ind)
490
       end do
491
       !$OMP END PARALLEL DO
492
#else
493
       T1V = alpha * T1V + omega * T2V
494
#endif
495
       residual = normfun( HUTI_NDIM, T1V, 1 )
496
    case (HUTI_USUPPLIED_STOPC)
497
       residual = stopcfun( X, B, R, ipar, dpar )
498
    case default
499
       call matvecsubr( X, T2V, ipar )
500
#ifdef _OPENMP
501
       !$OMP PARALLEL DO
502
       do i=1,ndim
503
          work(i,T1V_ind) = work(i,T2V_ind) - rhsvec(i)
504
       end do
505
       !$OMP END PARALLEL DO
506
#else
507
       T1V = T2V - B
508
#endif
509
       residual = normfun( HUTI_NDIM, T1V, 1 )
510
    end select
511

512
    !
513
    ! Print debugging output if desired
514
    !
515

516
    if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
517
       if ( mod(iter_count, HUTI_DBUGLVL) .eq. 0 ) then
518
          write (*, '(I8, E11.4)') iter_count, residual
519
       end if
520
    end if
521

522
    if ( residual .lt. HUTI_TOLERANCE ) then
523
       HUTI_INFO = HUTI_CONVERGENCE
524
       go to 1000
525
    end if
526

527
    if ( omega .eq. 0 ) then
528
       HUTI_INFO = HUTI_BICGSTAB_OMEGA
529
       go to 1000
530
    end if
531

532
    IF( residual /= residual .OR. residual > HUTI_MAXTOLERANCE ) THEN
533
      HUTI_INFO = HUTI_DIVERGENCE
534
      GOTO 1000
535
    END IF
536
    
537
    oldrho = rho
538

539
    !
540
    ! Return back to the iteration loop (without initialization)
541
    !
542

543
    iter_count = iter_count + 1
544
    if ( iter_count .gt. HUTI_MAXIT ) then
545
       HUTI_INFO = HUTI_MAXITER
546
       go to 1000
547
    end if
548

549
    go to 300
550

551
    !
552
    ! This is where we exit last time (after enough iterations or breakdown)
553
    !
554

555
1000 continue
556
    if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
557
       write (*, '(I8, E11.4)') iter_count, residual
558
    end if
559

560
    HUTI_ITERS = iter_count
561
    return
562

563
    ! End of execution
564
    !*********************************************************************
565

566
  end subroutine  huti_dbicgstabsolv
567

568
  !*************************************************************************
569

570

571
  !*************************************************************************
572
  !*************************************************************************
573
  ! Complex version
574
  !*************************************************************************
575
  !*************************************************************************
576

577
  subroutine  huti_cbicgstabsolv  ( ndim, wrkdim, xvec, rhsvec, &
578
       ipar, dpar, work, matvecsubr, pcondlsubr, &
579
       pcondrsubr, dotprodfun, normfun, stopcfun )
580

581
    use huti_interfaces
582
    implicit none
583

584
    procedure( mv_iface_c ), pointer :: matvecsubr
585
    procedure( pc_iface_c ), pointer :: pcondlsubr
586
    procedure( pc_iface_c ), pointer :: pcondrsubr
587
    procedure( dotp_iface_c ), pointer :: dotprodfun
588
    procedure( norm_iface_c ), pointer :: normfun
589
    procedure( stopc_iface_c ), pointer :: stopcfun
590

591
    ! Parameters
592

593
    integer :: ndim, wrkdim
594
    complex, dimension(ndim) :: xvec, rhsvec
595
    integer, dimension(HUTI_IPAR_DFLTSIZE) :: ipar
596
    double precision, dimension(HUTI_DPAR_DFLTSIZE) :: dpar
597
    complex, dimension(ndim,wrkdim) :: work
598

599
    ! Local variables
600

601
    complex :: rho, oldrho, alpha, beta, omega
602
    integer :: iter_count
603

604
    real :: residual, rhsnorm, precrhsnorm
605

606
    !
607
    ! End of variable declarations
608
    !*********************************************************************
609

610
    !*********************************************************************
611
    ! The actual BiCGSTAB begins here (look the pseudo code in the
612
    ! "Templates..."-book, page 27)
613
    !
614
    ! First the initialization part
615
    !
616

617
    iter_count = 1
618

619
    ! The following applies for all matrix operations in this solver
620

621
    HUTI_EXTOP_MATTYPE = HUTI_MAT_NOTTRPSED
622

623
    ! Norms of right-hand side vector are used in convergence tests
624

625
    if ( HUTI_STOPC .eq. HUTI_TRESID_SCALED_BYB .or. & 
626
         HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYB ) then
627
       rhsnorm = normfun( HUTI_NDIM, B, 1 )
628
    end if
629
    if ( HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYPRECB ) then
630
       call pcondlsubr( P, B, ipar )
631
       precrhsnorm = normfun( HUTI_NDIM, P, 1 )
632
    end if
633

634
    ! Generate vector X if needed
635

636
    if ( HUTI_INITIALX .eq. HUTI_RANDOMX ) then
637
       call  huti_crandvec   ( X, ipar )
638
    else if ( HUTI_INITIALX .ne. HUTI_USERSUPPLIEDX ) then
639
       X = 1
640
    end if
641

642
    call matvecsubr( X, R, ipar )
643
    R = B - R
644
    RTLD = R
645
    P = 0; V = 0
646
    oldrho = 1; omega = 1; alpha = 0
647

648
    !
649
    ! This is where the loop starts (that is we continue from here after
650
    ! the first iteration)
651
    !
652

653
300 continue
654

655
    rho = dotprodfun( HUTI_NDIM, RTLD, 1, R, 1 )
656
    if ( rho .eq. 0 ) then
657
       HUTI_INFO = HUTI_BICGSTAB_RHO
658
       go to 1000
659
    end if
660

661
    beta = ( rho * alpha ) / ( oldrho * omega )
662
    P = R + beta * ( P - omega * V )
663

664
    call pcondlsubr( V, P, ipar )
665
    call pcondrsubr( T1V, V, ipar )
666
    call matvecsubr( T1V, V, ipar )
667

668
    alpha = rho / dotprodfun( HUTI_NDIM, RTLD, 1, V, 1 )
669
    S = R - alpha * V
670

671
    residual = normfun( HUTI_NDIM, S, 1 )
672
    if ( residual .lt. HUTI_EPSILON ) then
673
       X = X + alpha * T1V
674
       HUTI_INFO = HUTI_CONVERGENCE
675
       !HUTI_INFO = HUTI_BICGSTAB_SNORM
676
       go to 1000
677
    end if
678

679
    call pcondlsubr( T, S, ipar )
680
    call pcondrsubr( T2V, T, ipar )
681
    call matvecsubr( T2V, T, ipar )
682

683
    omega = ( dotprodfun( HUTI_NDIM, T, 1, S, 1 ) ) / &
684
         ( dotprodfun( HUTI_NDIM, T, 1, T, 1 ) )
685
    X = X + alpha * T1V + omega * T2V
686
    R = S - omega * T
687

688
    !
689
    ! Check the convergence against selected stopping criterion
690
    !
691

692
    select case (HUTI_STOPC)
693
    case (HUTI_TRUERESIDUAL)
694
       call matvecsubr( X, T2V, ipar )
695
       T1V = T2V - B
696
       residual = normfun( HUTI_NDIM, T1V, 1 )
697
    case (HUTI_TRESID_SCALED_BYB)
698
       call matvecsubr( X, T2V, ipar )
699
       T1V = T2V - B
700
       residual = normfun( HUTI_NDIM, T1V, 1 ) / rhsnorm
701
    case (HUTI_PSEUDORESIDUAL)
702
       residual = normfun( HUTI_NDIM, R, 1 )
703
    case (HUTI_PRESID_SCALED_BYB)
704
       residual = normfun( HUTI_NDIM, R, 1 ) / rhsnorm
705
    case (HUTI_PRESID_SCALED_BYPRECB)
706
       residual = normfun( HUTI_NDIM, R, 1 ) / precrhsnorm
707
    case (HUTI_XDIFF_NORM)
708
       T1V = alpha * T1V + omega * T2V
709
       residual = normfun( HUTI_NDIM, T1V, 1 )
710
    case (HUTI_USUPPLIED_STOPC)
711
       residual = stopcfun( X, B, R, ipar, dpar )
712
    case default
713
       call matvecsubr( X, T2V, ipar )
714
       T1V = T2V - B
715
       residual = normfun( HUTI_NDIM, T1V, 1 )
716
    end select
717

718
    !
719
    ! Print debugging output if desired
720
    !
721

722
    if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
723
       if ( mod(iter_count, HUTI_DBUGLVL) .eq. 0 ) then
724
          write (*, '(I8, E11.4)') iter_count, residual
725
       end if
726
    end if
727

728
    if ( residual .lt. HUTI_TOLERANCE ) then
729
       HUTI_INFO = HUTI_CONVERGENCE
730
       go to 1000
731
    end if
732

733
    IF( residual /= residual .OR. residual > HUTI_MAXTOLERANCE ) THEN
734
      HUTI_INFO = HUTI_DIVERGENCE
735
      GOTO 1000
736
    END IF
737
    
738
    if ( omega .eq. 0 ) then
739
       HUTI_INFO = HUTI_BICGSTAB_OMEGA
740
       go to 1000
741
    end if
742

743
    oldrho = rho
744

745
    !
746
    ! Return back to the iteration loop (without initialization)
747
    !
748

749
    iter_count = iter_count + 1
750
    if ( iter_count .gt. HUTI_MAXIT ) then
751
       HUTI_INFO = HUTI_MAXITER
752
       go to 1000
753
    end if
754

755
    go to 300
756

757
    !
758
    ! This is where we exit last time (after enough iterations or breakdown)
759
    !
760

761
1000 continue
762
    if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
763
       write (*, '(I8, E11.4)') iter_count, residual
764
    end if
765

766
    HUTI_ITERS = iter_count
767
    return
768

769
    ! End of execution
770
    !*********************************************************************
771

772
  end subroutine  huti_cbicgstabsolv
773

774
  !*************************************************************************
775

776

777
  !*************************************************************************
778
  !*************************************************************************
779
  ! Double complex version
780
  !*************************************************************************
781
  !*************************************************************************
782

783
  subroutine  huti_zbicgstabsolv  ( ndim, wrkdim, xvec, rhsvec, &
784
       ipar, dpar, work, matvecsubr, pcondlsubr, &
785
       pcondrsubr, dotprodfun, normfun, stopcfun )
786

787
    use huti_interfaces
788
    implicit none
789

790
    procedure( mv_iface_z ), pointer :: matvecsubr
791
    procedure( pc_iface_z ), pointer :: pcondlsubr
792
    procedure( pc_iface_z ), pointer :: pcondrsubr
793
    procedure( dotp_iface_z ), pointer :: dotprodfun
794
    procedure( norm_iface_z ), pointer :: normfun
795
    procedure( stopc_iface_z ), pointer :: stopcfun
796

797
    ! Parameters
798

799
    integer :: ndim, wrkdim
800
    double complex, dimension(ndim) :: xvec, rhsvec
801
    integer, dimension(HUTI_IPAR_DFLTSIZE) :: ipar
802
    double precision, dimension(HUTI_DPAR_DFLTSIZE) :: dpar
803
    double complex, dimension(ndim,wrkdim) :: work
804

805
    ! Local variables
806

807
    double complex :: rho, oldrho, alpha, beta, omega
808
    integer :: iter_count
809

810
    double precision :: residual, rhsnorm, precrhsnorm
811

812
    !
813
    ! End of variable declarations
814
    !*********************************************************************
815

816
    !*********************************************************************
817
    ! The actual BiCGSTAB begins here (look the pseudo code in the
818
    ! "Templates..."-book, page 27)
819
    !
820
    ! First the initialization part
821
    !
822

823
    iter_count = 1
824

825
    ! The following applies for all matrix operations in this solver
826

827
    HUTI_EXTOP_MATTYPE = HUTI_MAT_NOTTRPSED
828

829
    ! Norms of right-hand side vector are used in convergence tests
830

831
    if ( HUTI_STOPC .eq. HUTI_TRESID_SCALED_BYB .or. & 
832
         HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYB ) then
833
       rhsnorm = normfun( HUTI_NDIM, B, 1 )
834
    end if
835
    if ( HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYPRECB ) then
836
       call pcondlsubr( P, B, ipar )
837
       precrhsnorm = normfun( HUTI_NDIM, P, 1 )
838
    end if
839

840
    ! Generate vector X if needed
841

842
    if ( HUTI_INITIALX .eq. HUTI_RANDOMX ) then
843
       call  huti_zrandvec   ( X, ipar )
844
    else if ( HUTI_INITIALX .ne. HUTI_USERSUPPLIEDX ) then
845
       X = 1
846
    end if
847

848
    call matvecsubr( X, R, ipar )
849
    R = B - R
850
    RTLD = R
851
    P = 0; V = 0
852
    oldrho = 1; omega = 1; alpha = 0
853

854
    !
855
    ! This is where the loop starts (that is we continue from here after
856
    ! the first iteration)
857
    !
858

859
300 continue
860

861
    rho = dotprodfun( HUTI_NDIM, RTLD, 1, R, 1 )
862
    if ( rho .eq. 0 ) then
863
       HUTI_INFO = HUTI_BICGSTAB_RHO
864
       go to 1000
865
    end if
866

867
    beta = ( rho * alpha ) / ( oldrho * omega )
868
    P = R + beta * ( P - omega * V )
869

870
    call pcondlsubr( V, P, ipar )
871
    call pcondrsubr( T1V, V, ipar )
872
    call matvecsubr( T1V, V, ipar )
873

874
    alpha = rho / dotprodfun( HUTI_NDIM, RTLD, 1, V, 1 )
875
    S = R - alpha * V
876

877
    residual = normfun( HUTI_NDIM, S, 1 )
878
    if ( residual .lt. HUTI_EPSILON ) then
879
       X = X + alpha * T1V
880
       !HUTI_INFO = HUTI_BICGSTAB_SNORM
881
       HUTI_INFO = HUTI_CONVERGENCE
882
       go to 1000
883
    end if
884

885
    call pcondlsubr( T, S, ipar )
886
    call pcondrsubr( T2V, T, ipar )
887
    call matvecsubr( T2V, T, ipar )
888

889
    omega = ( dotprodfun( HUTI_NDIM, T, 1, S, 1 ) ) / &
890
         ( dotprodfun( HUTI_NDIM, T, 1, T, 1 ) )
891
    X = X + alpha * T1V + omega * T2V
892
    R = S - omega * T
893

894
    !
895
    ! Check the convergence against selected stopping criterion
896
    !
897

898
    select case (HUTI_STOPC)
899
    case (HUTI_TRUERESIDUAL)
900
       call matvecsubr( X, T2V, ipar )
901
       T1V = T2V - B
902
       residual = normfun( HUTI_NDIM, T1V, 1 )
903
    case (HUTI_TRESID_SCALED_BYB)
904
       call matvecsubr( X, T2V, ipar )
905
       T1V = T2V - B
906
       residual = normfun( HUTI_NDIM, T1V, 1 ) / rhsnorm
907
    case (HUTI_PSEUDORESIDUAL)
908
       residual = normfun( HUTI_NDIM, R, 1 )
909
    case (HUTI_PRESID_SCALED_BYB)
910
       residual = normfun( HUTI_NDIM, R, 1 ) / rhsnorm
911
    case (HUTI_PRESID_SCALED_BYPRECB)
912
       residual = normfun( HUTI_NDIM, R, 1 ) / precrhsnorm
913
    case (HUTI_XDIFF_NORM)
914
       T1V = alpha * T1V + omega * T2V
915
       residual = normfun( HUTI_NDIM, T1V, 1 )
916
    case (HUTI_USUPPLIED_STOPC)
917
       residual = stopcfun( X, B, R, ipar, dpar )
918
    case default
919
       call matvecsubr( X, T2V, ipar )
920
       T1V = T2V - B
921
       residual = normfun( HUTI_NDIM, T1V, 1 )
922
    end select
923

924
    !
925
    ! Print debugging output if desired
926
    !
927

928
    if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
929
       if ( mod(iter_count, HUTI_DBUGLVL) .eq. 0 ) then
930
          write (*, '(I8, E11.4)') iter_count, residual
931
       end if
932
    end if
933

934
    if ( residual .lt. HUTI_TOLERANCE ) then
935
       HUTI_INFO = HUTI_CONVERGENCE
936
       go to 1000
937
    end if
938

939
    IF( residual /= residual .OR. residual > HUTI_MAXTOLERANCE ) THEN
940
      HUTI_INFO = HUTI_DIVERGENCE
941
      GOTO 1000
942
    END IF
943
    
944
    if ( omega .eq. 0 ) then
945
       HUTI_INFO = HUTI_BICGSTAB_OMEGA
946
       go to 1000
947
    end if
948

949
    oldrho = rho
950

951
    !
952
    ! Return back to the iteration loop (without initialization)
953
    !
954

955
    iter_count = iter_count + 1
956
    if ( iter_count .gt. HUTI_MAXIT ) then
957
       HUTI_INFO = HUTI_MAXITER
958
       go to 1000
959
    end if
960

961
    go to 300
962

963
    !
964
    ! This is where we exit last time (after enough iterations or breakdown)
965
    !
966

967
1000 continue
968
    if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
969
       write (*, '(I8, E11.4)') iter_count, residual
970
    end if
971

972
    HUTI_ITERS = iter_count
973
    return
974

975
    ! End of execution
976
    !*********************************************************************
977

978
  end subroutine  huti_zbicgstabsolv
979

980
  !*************************************************************************
981

982

983
end module huti_bicgstab
984

985
Product

Resources

Company
