CoCalc -- IterativeMethods.F90

GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/fem/src/IterativeMethods.F90
³²⁰³ views
1
!/*****************************************************************************/
2
! *
3
! *  Elmer, A Finite Element Software for Multiphysical Problems
4
! *
5
! *  Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland
6
! * 
7
! * This library is free software; you can redistribute it and/or
8
! * modify it under the terms of the GNU Lesser General Public
9
! * License as published by the Free Software Foundation; either
10
! * version 2.1 of the License, or (at your option) any later version.
11
! *
12
! * This library is distributed in the hope that it will be useful,
13
! * but WITHOUT ANY WARRANTY; without even the implied warranty of
14
! * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
15
! * Lesser General Public License for more details.
16
! * 
17
! * You should have received a copy of the GNU Lesser General Public
18
! * License along with this library (in file ../LGPL-2.1); if not, write 
19
! * to the Free Software Foundation, Inc., 51 Franklin Street, 
20
! * Fifth Floor, Boston, MA  02110-1301  USA
21
! *
22
! *****************************************************************************/
23
!
24
!/******************************************************************************
25
! *
26
! *  Authors: Juha Ruokolainen, Peter Råback, Mika Malinen, Martin van Gijzen
27
! *  Email:   [email protected]
28
! *  Web:     http://www.csc.fi/elmer
29
! *  Address: CSC - IT Center for Science Ltd.
30
! *           Keilaranta 14
31
! *           02101 Espoo, Finland 
32
! *
33
! *  Original Date: 20.9.2007
34
! *
35
! ****************************************************************************/
36

37
!> \ingroup ElmerLib 
38
!> \{
39

40
!------------------------------------------------------------------------------
41
!>  Module containing iterative methods. Uses the calling procedure of the 
42
!>  HUTIter package for similar interfacing. The idea is that the future 
43
!>  development of iterative methods could be placed in this module. 
44
!------------------------------------------------------------------------------
45

46

47
#include "huti_fdefs.h"
48

49
! if using old huti_fdefs.h, later obsolete
50
#ifndef HUTI_MAXTOLERANCE
51
#define HUTI_MAXTOLERANCE dpar(2)
52
#endif
53
#ifndef HUTI_SGSPARAM
54
#define HUTI_SGSPARAM dpar(3)
55
#endif
56
#ifndef HUTI_PSEUDOCOMPLEX
57
#define HUTI_PSEUDOCOMPLEX ipar(7)
58
#endif
59
#ifndef HUTI_BICGSTABL_L
60
#define HUTI_BICGSTABL_L ipar(16)
61
#endif
62
#ifndef HUTI_DIVERGENCE
63
#define HUTI_DIVERGENCE 3
64
#endif
65
#ifndef HUTI_GCR_RESTART
66
#define HUTI_GCR_RESTART ipar(17)
67
#endif
68
#ifndef HUTI_IDRS_S
69
#define HUTI_IDRS_S ipar(18)
70
#endif
71

72

73
MODULE IterativeMethods
74
  
75
  USE CRSMatrix  
76
  USE SParIterComm
77
  
78
  IMPLICIT NONE
79
  
80
  INTEGER :: nc
81
  LOGICAL :: Constrained
82

83
  TYPE(Matrix_t), POINTER, PRIVATE :: CM
84
  
85
CONTAINS
86
  
87

88
  ! When treating a complex system with iterative solver, norm and matrix-vector product are
89
  ! similar for real-valued and complex-valued systems. However, the inner product is different.
90
  ! For pseudo-complex systems this routine generates also the complex part of the product.
91
  ! This may have a favourable effect on convergence.
92
  !
93
  ! This routine has same API as the fully real-valued system but every second call returns
94
  ! the missing complex part.
95
  !
96
  ! This routine assumes that in x and y the values follow each other. 
97
  !-----------------------------------------------------------------------------------
98
  FUNCTION PseudoZDotProd( ndim, x, xind, y, yind ) RESULT( d )
99
  !-----------------------------------------------------------------------------------
100
    IMPLICIT NONE
101
    
102
    INTEGER :: ndim, xind, yind
103
    REAL(KIND=dp) :: x(*)
104
    REAL(KIND=dp) :: y(*)
105
    REAL(KIND=dp) :: d
106
        
107
    INTEGER :: i, callcount = 0
108
    REAL(KIND=dp) :: a,b
109
    
110
    SAVE callcount, a, b
111

112
    IF( callcount == 0 ) THEN    
113
      ! z = x^H*y = (x_re-i*x_im)(y_re+i*y_im)       
114
      ! =>  z_re = x_re*y_re + x_im*y_im
115
      !     z_im = x_re*y_im - x_im*y_re
116
      
117
      a = SUM( x(1:ndim) * y(1:ndim) )
118
      b = SUM( x(1:ndim:2) * y(2:ndim:2) - x(2:ndim:2) * y(1:ndim:2) )
119

120
      IF (ParEnv % PEs > 1) THEN
121
        CALL SParActiveSUM(a,0)
122
        CALL SParActiveSUM(b,0)
123
      END IF
124

125
      d = a 
126
      callcount = callcount + 1
127
    ELSE
128
      d = b
129
      callcount = 0
130
    END IF
131
      
132
    !-----------------------------------------------------------------------------------
133
  END FUNCTION PseudoZDotProd
134
  !-----------------------------------------------------------------------------------
135

136
  
137
  ! As the previous but assumes that the real and complex values are ordered blockwise. 
138
  !-----------------------------------------------------------------------------------
139
  FUNCTION PseudoZDotProd2( ndim, x, xind, y, yind ) RESULT( d )
140
  !-----------------------------------------------------------------------------------
141
    IMPLICIT NONE
142
    
143
    INTEGER :: ndim, xind, yind
144
    REAL(KIND=dp) :: x(*)
145
    REAL(KIND=dp) :: y(*)
146
    REAL(KIND=dp) :: d
147
        
148
    INTEGER :: i, callcount = 0
149
    REAL(KIND=dp) :: a,b
150
    
151
    SAVE callcount, a, b
152

153
    IF( callcount == 0 ) THEN    
154
      a = SUM( x(1:ndim) * y(1:ndim) )
155
      b = SUM( x(1:ndim/2) * y(ndim/2+1:ndim) - x(ndim/2+1:ndim) * y(1:ndim/2) )
156

157
      IF (ParEnv % PEs > 1) THEN
158
        CALL SParActiveSUM(a,0)
159
        CALL SParActiveSUM(b,0)
160
      END IF
161
      
162
      d = a 
163
      callcount = callcount + 1
164
    ELSE
165
      d = b
166
      callcount = 0
167
    END IF
168
      
169
!-----------------------------------------------------------------------------------
170
  END FUNCTION PseudoZDotProd2
171
!-----------------------------------------------------------------------------------
172

173
  
174
!------------------------------------------------------------------------------
175
!> Symmetric Gauss-Seidel iterative method for linear systems. This is not really of practical
176
!> use but may be used for testing, for example. 
177
!------------------------------------------------------------------------------
178
  SUBROUTINE itermethod_sgs( xvec, rhsvec, &
179
      ipar, dpar, work, matvecsubr, pcondlsubr, &
180
      pcondrsubr, dotprodfun, normfun, stopcfun )
181

182
    USE huti_interfaces
183
    IMPLICIT NONE
184
    PROCEDURE( mv_iface_d ), POINTER :: matvecsubr
185
    PROCEDURE( pc_iface_d ), POINTER :: pcondlsubr
186
    PROCEDURE( pc_iface_d ), POINTER :: pcondrsubr
187
    PROCEDURE( dotp_iface_d ), POINTER :: dotprodfun
188
    PROCEDURE( norm_iface_d ), POINTER :: normfun
189
    PROCEDURE( stopc_iface_d ), POINTER :: stopcfun
190

191
    INTEGER, DIMENSION(HUTI_IPAR_DFLTSIZE) :: ipar
192
    REAL(KIND=dp), DIMENSION(HUTI_NDIM) :: xvec, rhsvec
193
    DOUBLE PRECISION, DIMENSION(HUTI_DPAR_DFLTSIZE) :: dpar
194
    REAL(KIND=dp), DIMENSION(HUTI_WRKDIM,HUTI_NDIM) :: work
195

196
    INTEGER :: ndim
197
    INTEGER :: Rounds, OutputInterval
198
    REAL(KIND=dp) :: MinTol, MaxTol, Residual, Omega
199
    LOGICAL :: Converged, Diverged
200

201
    ndim = HUTI_NDIM
202
    Rounds = HUTI_MAXIT
203
    MinTol = HUTI_TOLERANCE
204
    MaxTol = HUTI_MAXTOLERANCE
205
    OutputInterval = HUTI_DBUGLVL 
206
    Omega = HUTI_SGSPARAM
207

208
    CALL sgs(ndim, GlobalMatrix, xvec, rhsvec, Rounds, MinTol, MaxTol, Residual, &
209
        Converged, Diverged, OutputInterval, Omega)
210

211
    IF(Converged) HUTI_INFO = HUTI_CONVERGENCE
212
    IF(Diverged) HUTI_INFO = HUTI_DIVERGENCE
213
    IF ( (.NOT. Converged) .AND. (.NOT. Diverged) ) HUTI_INFO = HUTI_MAXITER
214
    
215
  CONTAINS 
216
 
217
!------------------------------------------------------------------------------
218
    SUBROUTINE SGS( n, A, x, b, Rounds, MinTolerance, MaxTolerance, Residual, &
219
        Converged, Diverged, OutputInterval, Omega )
220
 !------------------------------------------------------------------------------
221
      TYPE(Matrix_t), POINTER :: A
222
      INTEGER :: Rounds
223
      INTEGER :: i,j,k,n
224
      REAL(KIND=dp) :: x(n),b(n)
225
      REAL(KIND=dp) :: Omega
226
      INTEGER, POINTER :: Cols(:),Rows(:)
227
      REAL(KIND=dp), POINTER :: Values(:)
228
      REAL(KIND=dp), ALLOCATABLE :: R(:)
229
      LOGICAL :: Converged, Diverged
230
      INTEGER :: OutputInterval
231
      REAL(KIND=dp) :: MinTolerance, MaxTolerance, Residual, bnorm,rnorm,w,s
232

233
      Rows   => A % Rows
234
      Cols   => A % Cols
235
      Values => A % Values
236
      
237
      ALLOCATE( R(n) )
238
      
239
      CALL matvecsubr( x, r, ipar )
240
     
241
      r(1:n) = b(1:n) - r(1:n)
242
      bnorm = normfun(n, b, 1)
243
      rnorm = normfun(n, r, 1) 
244

245
      Residual = rnorm / bnorm
246
      Converged = (Residual < MinTolerance) 
247
      Diverged = (Residual > MaxTolerance) .OR. (Residual /= Residual)
248
      IF( Converged .OR. Diverged) RETURN
249

250
      DO k=1,Rounds
251
        DO i=1,n
252
          s = 0.0d0
253
          DO j=Rows(i),Rows(i+1)-1
254
            s = s + x(Cols(j)) * Values(j)
255
          END DO
256
          x(i) = x(i) + Omega * (b(i)-s) / Values(A % Diag(i))
257
        END DO
258
        
259
        DO i=n,1,-1
260
          s = 0.0d0
261
          DO j=Rows(i),Rows(i+1)-1
262
            s = s + x(Cols(j)) * Values(j)
263
          END DO
264
          x(i) = x(i) + Omega * (b(i)-s) / Values(A % Diag(i))
265
        END DO
266
        
267
        CALL matvecsubr( x, r, ipar )
268
        r(1:n) = b(1:n) - r(1:n)
269
        rnorm = normfun(n, r, 1)
270
        
271
        Residual = rnorm / bnorm
272
        IF( MOD(k,OutputInterval) == 0) THEN
273
          WRITE (*, '(I8, 2E11.4)') k, rnorm, residual
274
          CALL FLUSH(6)
275
        END IF
276
        
277
        Converged = (Residual < MinTolerance) 
278
        Diverged = (Residual > MaxTolerance) .OR. (Residual /= Residual)
279
        IF( Converged .OR. Diverged) RETURN
280
        
281
      END DO
282
    END SUBROUTINE SGS
283
!------------------------------------------------------------------------------
284
  END SUBROUTINE itermethod_sgs
285
 !------------------------------------------------------------------------------
286
 
287

288

289
!------------------------------------------------------------------------------
290
!> Jacobi iterative method for linear systems. This is not really of practical
291
!> use but may be used for testing, for example. 
292
!> Note that if the scaling is performed so that the diagonal entry is one
293
!> the division by it is unnecessary. Hence for this method scaling is not 
294
!> needed. 
295
!------------------------------------------------------------------------------
296
 SUBROUTINE itermethod_jacobi( xvec, rhsvec, &
297
      ipar, dpar, work, matvecsubr, pcondlsubr, &
298
      pcondrsubr, dotprodfun, normfun, stopcfun )
299

300
    USE huti_interfaces
301
    IMPLICIT NONE
302
    PROCEDURE( mv_iface_d ), POINTER :: matvecsubr
303
    PROCEDURE( pc_iface_d ), POINTER :: pcondlsubr
304
    PROCEDURE( pc_iface_d ), POINTER :: pcondrsubr
305
    PROCEDURE( dotp_iface_d ), POINTER :: dotprodfun
306
    PROCEDURE( norm_iface_d ), POINTER :: normfun
307
    PROCEDURE( stopc_iface_d ), POINTER :: stopcfun
308

309
    INTEGER, DIMENSION(HUTI_IPAR_DFLTSIZE) :: ipar
310
    REAL(KIND=dp), TARGET, DIMENSION(HUTI_NDIM) :: xvec, rhsvec
311
    DOUBLE PRECISION, DIMENSION(HUTI_DPAR_DFLTSIZE) :: dpar
312
    REAL(KIND=dp), DIMENSION(HUTI_WRKDIM,HUTI_NDIM) :: work
313

314
    INTEGER :: ndim
315
    INTEGER :: Rounds, OutputInterval
316
    REAL(KIND=dp) :: MinTol, MaxTol, Residual
317
    LOGICAL :: Converged, Diverged
318
    
319
    ndim = HUTI_NDIM
320
    Rounds = HUTI_MAXIT
321
    MinTol = HUTI_TOLERANCE
322
    MaxTol = HUTI_MAXTOLERANCE
323
    OutputInterval = HUTI_DBUGLVL 
324
       
325
    CALL jacobi(ndim, GlobalMatrix, xvec, rhsvec, Rounds, MinTol, MaxTol, Residual, &
326
        Converged, Diverged, OutputInterval )
327

328
    IF(Converged) HUTI_INFO = HUTI_CONVERGENCE
329
    IF(Diverged) HUTI_INFO = HUTI_DIVERGENCE
330
    IF ( (.NOT. Converged) .AND. (.NOT. Diverged) ) HUTI_INFO = HUTI_MAXITER   
331

332
  CONTAINS 
333
    
334
    
335
    SUBROUTINE Jacobi( n, A, x, b, Rounds, MinTolerance, MaxTolerance, Residual, &
336
        Converged, Diverged, OutputInterval) 
337
!------------------------------------------------------------------------------
338
      TYPE(Matrix_t), POINTER :: A
339
      INTEGER :: Rounds
340
      REAL(KIND=dp) :: x(n),b(n)
341
      LOGICAL :: Converged, Diverged
342
      REAL(KIND=dp) :: MinTolerance, MaxTolerance, Residual
343
      INTEGER :: OutputInterval
344
      REAL(KIND=dp) :: bnorm,rnorm
345
      REAL(KIND=dp), ALLOCATABLE :: R(:)
346
!------------------------------------------------------------------------------
347
      INTEGER :: i,j,n
348
!------------------------------------------------------------------------------
349
      
350
      Converged = .FALSE.
351
      Diverged = .FALSE.
352
      
353
      ALLOCATE( R(n) )
354
      
355
      CALL matvecsubr( x, r, ipar )
356
      r(1:n) = b(1:n) - r(1:n)
357
      
358
      bnorm = normfun(n, b, 1)
359
      rnorm = normfun(n, r, 1)
360
      
361
      Residual = rnorm / bnorm
362
      Converged = (Residual < MinTolerance) 
363
      Diverged = (Residual > MaxTolerance) .OR. (Residual /= Residual)    
364
      IF( Converged .OR. Diverged) RETURN
365
      
366
      DO i=1,Rounds
367
        DO j=1,n
368
          x(j) = x(j) + r(j) / A % Values(A % diag(j))
369
        END DO
370
        CALL matvecsubr( x, r, ipar )
371
        
372
        r(1:n) = b(1:n) - r(1:n)
373
        rnorm = normfun(n, r, 1)
374
        
375
        Residual = rnorm / bnorm
376
        
377
        IF( MOD(i,OutputInterval) == 0) THEN
378
          WRITE (*, '(I8, 2E11.4)') i, rnorm, residual
379
          CALL FLUSH(6)
380
        END IF
381
        
382
        Converged = (Residual < MinTolerance) 
383
        Diverged = (Residual > MaxTolerance) .OR. (Residual /= Residual)
384
        IF( Converged .OR. Diverged) EXIT
385
      END DO
386
      
387
      DEALLOCATE( R )
388
      
389
    END SUBROUTINE Jacobi
390
    
391
!------------------------------------------------------------------------------
392
  END SUBROUTINE itermethod_jacobi
393
!------------------------------------------------------------------------------
394
  
395

396
!------------------------------------------------------------------------------
397
!> Richardson iterative method for linear systems. This may of actual use for 
398
!> mass matrices. Actually this is not the simple Richardson iteration method
399
!> as it is preconditioned with the lumped mass matrix. 
400
!> Note that if scaling is performed by the "row equilibrium" method then
401
!> lumped mass is by construction unity (assuming all-positive entries).
402
!> So for this method scaling is not needed. 
403
!------------------------------------------------------------------------------
404
 SUBROUTINE itermethod_richardson( xvec, rhsvec, &
405
      ipar, dpar, work, matvecsubr, pcondlsubr, &
406
      pcondrsubr, dotprodfun, normfun, stopcfun )
407

408
    USE huti_interfaces
409
    IMPLICIT NONE
410
    PROCEDURE( mv_iface_d ), POINTER :: matvecsubr
411
    PROCEDURE( pc_iface_d ), POINTER :: pcondlsubr
412
    PROCEDURE( pc_iface_d ), POINTER :: pcondrsubr
413
    PROCEDURE( dotp_iface_d ), POINTER :: dotprodfun
414
    PROCEDURE( norm_iface_d ), POINTER :: normfun
415
    PROCEDURE( stopc_iface_d ), POINTER :: stopcfun
416

417
    INTEGER, DIMENSION(HUTI_IPAR_DFLTSIZE) :: ipar
418
    REAL(KIND=dp), TARGET, DIMENSION(HUTI_NDIM) :: xvec, rhsvec
419
    DOUBLE PRECISION, DIMENSION(HUTI_DPAR_DFLTSIZE) :: dpar
420
    REAL(KIND=dp), DIMENSION(HUTI_WRKDIM,HUTI_NDIM) :: work
421

422
    INTEGER :: ndim
423
    INTEGER :: Rounds, OutputInterval
424
    REAL(KIND=dp) :: MinTol, MaxTol, Residual
425
    LOGICAL :: Converged, Diverged
426
    
427
    ndim = HUTI_NDIM
428
    Rounds = HUTI_MAXIT
429
    MinTol = HUTI_TOLERANCE
430
    MaxTol = HUTI_MAXTOLERANCE
431
    OutputInterval = HUTI_DBUGLVL 
432
       
433
    CALL richardson(ndim, GlobalMatrix, xvec, rhsvec, Rounds, MinTol, MaxTol, Residual, &
434
        Converged, Diverged, OutputInterval )
435

436
    IF(Converged) HUTI_INFO = HUTI_CONVERGENCE
437
    IF(Diverged) HUTI_INFO = HUTI_DIVERGENCE
438
    IF ( (.NOT. Converged) .AND. (.NOT. Diverged) ) HUTI_INFO = HUTI_MAXITER   
439

440
  CONTAINS 
441
    
442
    
443
    SUBROUTINE Richardson( n, A, x, b, Rounds, MinTolerance, MaxTolerance, Residual, &
444
        Converged, Diverged, OutputInterval) 
445
!------------------------------------------------------------------------------
446
      TYPE(Matrix_t), POINTER :: A
447
      INTEGER :: Rounds
448
      REAL(KIND=dp) :: x(n),b(n)
449
      LOGICAL :: Converged, Diverged
450
      REAL(KIND=dp) :: MinTolerance, MaxTolerance, Residual
451
      INTEGER :: OutputInterval
452
      REAL(KIND=dp) :: bnorm,rnorm,s,q
453
      INTEGER, POINTER :: Cols(:),Rows(:)
454
      REAL(KIND=dp), POINTER :: Values(:)
455
      REAL(KIND=dp), ALLOCATABLE :: R(:), M(:)
456
!------------------------------------------------------------------------------
457
      INTEGER :: i,j,k,n
458
!------------------------------------------------------------------------------
459

460
      Rows   => A % Rows
461
      Cols   => A % Cols
462
      Values => A % Values
463
      
464
      Converged = .FALSE.
465
      Diverged = .FALSE.
466
      
467
      ALLOCATE( R(n), M(n) )
468
      
469
      CALL matvecsubr( x, r, ipar )
470
      r(1:n) = b(1:n) - r(1:n)
471
      
472
      bnorm = normfun(n, b, 1)
473
      rnorm = normfun(n, r, 1)
474

475
      Residual = rnorm / bnorm
476
      Converged = (Residual < MinTolerance) 
477
      Diverged = (Residual > MaxTolerance) .OR. (Residual /= Residual)
478
      IF( Converged .OR. Diverged) RETURN
479

480
      ! Perform preconditioning by mass lumping
481
      DO i=1,n
482
        s = 0.0_dp
483
        DO j=Rows(i),Rows(i+1)-1
484
          s = s + Values( j )
485
        END DO
486
        M(i) = s 
487
      END DO
488

489
      DO k=1,Rounds
490
        DO i=1,n
491
          IF( k == 1 ) THEN
492
            x(i) = b(i) / M(i) 
493
          ELSE
494
            x(i) = x(i) + r(i) / M(i)
495
          END IF
496
        END DO
497
        
498
        CALL matvecsubr( x, r, ipar )
499

500
        r(1:n) = b(1:n) - r(1:n)
501
        rnorm = normfun(n, r, 1)
502
        
503
        Residual = rnorm / bnorm
504
        
505
        IF( MOD(k,OutputInterval) == 0) THEN
506
          WRITE (*, '(I8, 2E11.4)') k, rnorm, residual
507
          CALL FLUSH(6)
508
        END IF
509
        
510
        Converged = (Residual < MinTolerance) 
511
        Diverged = (Residual > MaxTolerance) .OR. (Residual /= Residual)
512
        IF( Converged .OR. Diverged) EXIT
513
      END DO
514
      
515
      DEALLOCATE( R, M )
516
      
517
    END SUBROUTINE Richardson
518
    
519
!------------------------------------------------------------------------------
520
  END SUBROUTINE itermethod_richardson
521
!------------------------------------------------------------------------------
522
  
523

524
!-----------------------------------------------------------------------------------
525
    SUBROUTINE C_matvec(u,v,ipar,matvecsubr)
526
!-----------------------------------------------------------------------------------
527
      USE huti_interfaces
528
      IMPLICIT NONE
529
      PROCEDURE( mv_iface_d ), POINTER :: matvecsubr
530
      INTEGER :: ipar(*)
531
      REAL(KIND=dp) :: u(*),v(*)
532
!-----------------------------------------------------------------------------------
533
      INTEGER :: i,j,k,l,ndim
534
!-----------------------------------------------------------------------------------
535
      ndim = HUTI_NDIM
536
      CALL matvecsubr( u, v, ipar )
537
      IF (Constrained) THEN
538
         DO i=1,CM % NumberOfRows
539
           k = ndim+i
540
           v(k) = 0._dp
541
           DO j = CM % Rows(i),CM % Rows(i+1)-1
542
             l = CM % Cols(j)
543
             v(l) = v(l) + CM % Values(j)*u(k)
544
             v(k) = v(k) + CM % Values(j)*u(l)
545
           END DO
546
         END DO
547
       END IF
548
!-----------------------------------------------------------------------------------
549
    END SUBROUTINE C_matvec
550
!-----------------------------------------------------------------------------------
551

552

553
!-----------------------------------------------------------------------------------
554
    RECURSIVE SUBROUTINE C_rpcond(u,v,ipar,pcondrsubr)
555
!-----------------------------------------------------------------------------------
556
      USE huti_interfaces
557
      IMPLICIT NONE
558
      PROCEDURE( pc_iface_d ), POINTER :: pcondrsubr
559
      INTEGER :: ipar(*)
560
      REAL(KIND=dp) :: u(*),v(*)
561

562
!-----------------------------------------------------------------------------------
563
      INTEGER :: ndim
564
!-----------------------------------------------------------------------------------
565
      ndim = HUTI_NDIM
566
      IF(Constrained) HUTI_NDIM = ndim+nc
567
      CALL pcondrsubr(u,v,ipar)
568
      IF(Constrained) HUTI_NDIM = ndim
569
!-----------------------------------------------------------------------------------
570
    END SUBROUTINE C_rpcond
571
!-----------------------------------------------------------------------------------
572

573
!------------------------------------------------------------------------------
574
!>   This routine solves real linear systems Ax = b by using the BiCGStab(l) algorithm 
575
!>   with l >= 2 and the right-oriented ILU(n) preconditioning. 
576
!------------------------------------------------------------------------------
577
  SUBROUTINE itermethod_bicgstabl( xvec, rhsvec, &
578
      ipar, dpar, work, matvecsubr, pcondlsubr, &
579
      pcondrsubr, dotprodfun, normfun, stopcfun )
580

581
    USE huti_interfaces
582
    IMPLICIT NONE
583
    PROCEDURE( mv_iface_d ), POINTER :: matvecsubr
584
    PROCEDURE( pc_iface_d ), POINTER :: pcondlsubr
585
    PROCEDURE( pc_iface_d ), POINTER :: pcondrsubr
586
    PROCEDURE( dotp_iface_d ), POINTER :: dotprodfun
587
    PROCEDURE( norm_iface_d ), POINTER :: normfun
588
    PROCEDURE( stopc_iface_d ), POINTER :: stopcfun
589

590
    INTEGER, DIMENSION(HUTI_IPAR_DFLTSIZE) :: ipar
591
    REAL(KIND=dp), TARGET, DIMENSION(HUTI_NDIM) :: xvec, rhsvec
592
    ! DOUBLE PRECISION, DIMENSION(HUTI_NDIM), TARGET :: xvec, rhsvec
593
    DOUBLE PRECISION, DIMENSION(HUTI_DPAR_DFLTSIZE) :: dpar
594
    REAL(KIND=dp), DIMENSION(HUTI_WRKDIM,HUTI_NDIM) :: work
595
    ! DOUBLE PRECISION, DIMENSION(HUTI_WRKDIM,HUTI_NDIM) :: work
596

597
    INTEGER :: ndim,i,j,k
598
    INTEGER :: Rounds, OutputInterval, PolynomialDegree
599
    REAL(KIND=dp) :: MinTol, MaxTol, Residual
600
    LOGICAL :: Converged, Diverged, Halted, UseStopCFun, PseudoComplex
601

602
    TYPE(Matrix_t),POINTER :: A
603

604
    REAL(KIND=dp), POINTER CONTIG :: x(:),b(:)
605

606
    ! Variables related to robust mode
607
    LOGICAL :: Robust 
608
    INTEGER :: BestIter,BadIterCount,MaxBadIter, RobustStart
609
    REAL(KIND=dp) :: BestNorm,RobustStep,RobustTol,RobustMaxTol
610
    REAL(KIND=dp), ALLOCATABLE :: Bestx(:)
611

612
    
613
    A => GlobalMatrix
614
    CM => A % ConstraintMatrix
615
    Constrained = ASSOCIATED(CM)
616
    
617
    ndim = HUTI_NDIM
618

619
    x => xvec
620
    b => rhsvec
621
    nc = 0
622
    IF (Constrained) THEN
623
      nc = CM % NumberOfRows
624
      Constrained = nc>0
625
      IF(Constrained) THEN
626
        ALLOCATE(x(ndim+nc),b(ndim+nc))
627
        IF(.NOT.ALLOCATED(CM % ExtraVals))THEN
628
          ALLOCATE(CM % ExtraVals(nc)); CM % extraVals=0._dp
629
        END IF
630
        b(1:ndim) = rhsvec; b(ndim+1:) = CM % RHS
631
        x(1:ndim) = xvec; x(ndim+1:) = CM % extraVals
632
      END IF
633
    END IF
634

635
    Rounds = HUTI_MAXIT
636
    MinTol = HUTI_TOLERANCE
637
    MaxTol = HUTI_MAXTOLERANCE
638
    OutputInterval = HUTI_DBUGLVL 
639
    PolynomialDegree = HUTI_BICGSTABL_L 
640
    UseStopCFun = HUTI_STOPC == HUTI_USUPPLIED_STOPC
641

642
    PseudoComplex = ( HUTI_PSEUDOCOMPLEX > 0 )  
643
    
644
    Converged = .FALSE.
645
    Diverged = .FALSE.
646
    Halted = .FALSE.
647
    
648
    Robust = ( HUTI_ROBUST == 1 )
649
    IF( Robust ) THEN
650
      RobustTol = HUTI_ROBUST_TOLERANCE
651
      RobustStep = HUTI_ROBUST_STEPSIZE
652
      RobustMaxTol = HUTI_ROBUST_MAXTOLERANCE
653
      MaxBadIter = HUTI_ROBUST_MAXBADIT
654
      RobustStart = HUTI_ROBUST_START
655
      BestNorm = SQRT(HUGE(BestNorm))
656
      BadIterCount = 0
657
      BestIter = 0      
658
      ALLOCATE( BestX(ndim))
659
    END IF
660
    
661
    CALL RealBiCGStabl(ndim+nc, A,x,b, Rounds, MinTol, MaxTol, &
662
         Converged, Diverged, Halted, OutputInterval, PolynomialDegree )
663

664
    IF(Constrained) THEN
665
      xvec=x(1:ndim)
666
      rhsvec=b(1:ndim)
667
      CM % extraVals = x(ndim+1:ndim+nc)
668
      DEALLOCATE(x,b)
669
    END IF
670
    
671
    IF( Robust ) THEN
672
      DEALLOCATE( BestX )
673
    END IF
674
      
675
    IF(Converged) THEN
676
      HUTI_INFO = HUTI_CONVERGENCE
677
    ELSE IF(Diverged) THEN
678
      HUTI_INFO = HUTI_DIVERGENCE
679
    ELSE IF(Halted) THEN
680
      HUTI_INFO = HUTI_HALTED
681
    ELSE
682
      HUTI_INFO = HUTI_MAXITER
683
    END IF
684
      
685
  CONTAINS
686

687
!-----------------------------------------------------------------------------------
688
!>   The subroutine has been written using as a starting point the work of D.R. Fokkema 
689
!>   (subroutine zbistbl v1.1 1998). Dr. Fokkema has given the right to distribute
690
!>   the derived work under GPL and hence the original more conservative 
691
!> copyright notice of the subroutine has been removed accordingly.  
692
!-----------------------------------------------------------------------------------
693
    SUBROUTINE RealBiCGStabl( n, A, x, b, MaxRounds, Tol, MaxTol, Converged, &
694
        Diverged, Halted, OutputInterval, l)
695
!----------------------------------------------------------------------------------- 
696
      INTEGER :: l   ! polynomial degree
697
      INTEGER :: n, MaxRounds, OutputInterval   
698
      LOGICAL :: Converged, Diverged, Halted
699
      TYPE(Matrix_t), POINTER :: A
700
      REAL(KIND=dp) :: x(n), b(n)
701
      REAL(KIND=dp) :: Tol, MaxTol
702
!------------------------------------------------------------------------------
703
      REAL(KIND=dp) :: zero, one, t(n), kappa0, kappal 
704
      REAL(KIND=dp) :: dnrm2, rnrm0, rnrm, mxnrmx, mxnrmr, errorind, &
705
          delta = 1.0d-2, bnrm, bw_errorind, tottime
706
      INTEGER :: i, j, rr, r, u, xp, bp, z, zz, y0, yl, y, k, iwork(l-1), stat, Round, &
707
          IluOrder
708
      REAL(KIND=dp) :: alpha, beta, omega, rho0, rho1, sigma, ddot, varrho, hatgamma
709
      LOGICAL rcmp, xpdt, GotIt, EarlyExit
710
      REAL(KIND=dp), ALLOCATABLE :: work(:,:)
711
      REAL(KIND=dp) :: rwork(l+1,3+2*(l+1))
712
      REAL(KIND=dp) :: tmpmtr(l-1,l-1), tmpvec(l-1)
713
      REAL(KIND=dp) :: beta_im
714
!------------------------------------------------------------------------------
715
    
716
      IF ( l < 2) CALL Fatal( 'RealBiCGStabl', 'Polynomial degree < 2' )
717
      
718
      IF ( ALL(x == 0.0d0) ) x = b
719

720
      zero = 0.0d0
721
      one =  1.0d0
722

723
      ALLOCATE( work(n,3+2*(l+1)) )
724
      !$OMP PARALLEL PRIVATE(j)
725
      DO j=1,3+2*(l+1)
726
         !$OMP DO SCHEDULE(STATIC)
727
         DO i=1,n
728
            work(i,j) = 0.0d0
729
         END DO
730
         !$OMP END DO
731
      END DO
732
      !$OMP END PARALLEL
733
      DO j=1,3+2*(l+1)
734
         DO i=1,l+1
735
            rwork(i,j) = 0.0d0
736
         END DO
737
      END DO
738

739
      rr = 1
740
      r = rr+1
741
      u = r+(l+1)
742
      xp = u+(l+1)
743
      bp = xp+1
744
    
745
      z = 1
746
      zz = z+(l+1)
747
      y0 = zz+(l+1)
748
      yl = y0+1
749
      y = yl+1
750
    
751
      ! CALL C_matvec(x,work(:,r),ipar,matvecsubr)
752
      CALL C_matvec(x,work(1,r),ipar,matvecsubr)
753
      
754
      !$OMP PARALLEL DO SCHEDULE(STATIC)
755
      DO i=1,n
756
         work(i,r) = b(i) - work(i,r)
757
      END DO
758
      !$OMP END PARALLEL DO
759
      ! bnrm  = normfun(n, b(1:n), 1 )
760
      ! rnrm0 = normfun(n, work(1:n,r), 1 )
761
      bnrm  = normfun(n, b(1), 1 )
762
      rnrm0 = normfun(n, work(1,r), 1 )
763

764
      !-------------------------------------------------------------------
765
      ! Check whether the initial guess is already converged, diverged or NaN
766
      !--------------------------------------------------------------------
767
      Diverged = .FALSE.
768
      IF (bnrm /= bnrm) THEN
769
        CALL Fatal( 'RealBiCGStab(l)', 'Breakdown error: bnrm = NaN.' )
770
      ENDIF
771
      IF(rnrm0 /= rnrm0 ) THEN
772
        CALL Fatal( 'RealBiCGStab(l)', 'Breakdown error: nrm0 = NaN.' )
773
      END IF
774

775
      errorind = rnrm0 / bnrm
776
      IF(errorind /= errorind ) THEN
777
        CALL Fatal( 'RealBiCGStab(l)', 'Breakdown error: errorind = NaN.' )
778
      END IF
779
     
780
      Converged = (errorind < Tol)
781
      Diverged = (errorind > MaxTol) 
782

783
      IF( Converged .OR. Diverged ) RETURN
784

785
      EarlyExit = .FALSE.
786

787
      !$OMP PARALLEL 
788
      !$OMP DO SCHEDULE(STATIC)
789
      DO i=1,n
790
         work(i,rr) = work(i,r) 
791
         work(i,bp) = work(i,r)
792
      END DO
793
      !$OMP END DO NOWAIT
794
      !$OMP DO SCHEDULE(STATIC)
795
      DO i=1,n
796
         work(i,xp) = x(i)
797
      END DO
798
      !$OMP END DO NOWAIT
799
      !$OMP DO SCHEDULE(STATIC)
800
      DO i=1,n
801
         x(i) = zero
802
      END DO
803
      !$OMP END DO NOWAIT
804
      !$OMP END PARALLEL
805

806
      rnrm = rnrm0
807
      mxnrmx = rnrm0
808
      mxnrmr = rnrm0  
809
      alpha = zero
810
      omega = one
811
      sigma = one
812
      rho0 = one
813

814
      DO Round=1,MaxRounds
815
        !-------------------------
816
        ! --- The BiCG part ---
817
        !-------------------------
818
        rho0 = -omega*rho0
819
        
820
        DO k=1,l
821
          ! rho1 = dotprodfun(n, work(1:n,rr), 1, work(1:n,r+k-1), 1 )
822
          rho1 = dotprodfun(n, work(1,rr), 1, work(1,r+k-1), 1 )
823
          IF (rho0 == zero) THEN
824
            CALL Warn( 'RealBiCGStab(l)', 'Iteration halted: rho0 == zero.' )
825
            Halted = .TRUE.
826
            GOTO 100
827
          ENDIF
828
          IF (rho1 /= rho1) THEN
829
            CALL Fatal( 'RealBiCGStab(l)', 'Breakdown error: rho1 == NaN.' )
830
          ENDIF
831
         
832
          beta = alpha*(rho1/rho0)
833
          rho0 = rho1
834
          !$OMP PARALLEL PRIVATE(j)
835
          DO j=0,k-1
836
             !$OMP DO SCHEDULE(STATIC)
837
             DO i=1,n
838
                work(i,u+j) = work(i,r+j) - beta*work(i,u+j)
839
             END DO
840
             !$OMP END DO
841
          ENDDO
842
          !$OMP END PARALLEL
843

844
          ! CALL C_rpcond( t, work(:,u+k-1), ipar, pcondrsubr )
845
          CALL C_rpcond( t, work(1,u+k-1), ipar, pcondrsubr )
846
          ! CALL C_matvec( t, work(:,u+k), ipar, matvecsubr )
847
          CALL C_matvec( t, work(1,u+k), ipar, matvecsubr )
848
          ! sigma = dotprodfun(n, work(1:n,rr), 1, work(1:n,u+k), 1 )
849
          sigma = dotprodfun(n, work(1,rr), 1, work(1,u+k), 1 )
850
          
851
          IF (sigma == zero) THEN
852
            CALL Warn( 'RealBiCGStab(l)', 'Iteration halted: sigma == zero.' )
853
            Halted = .TRUE.
854
            GOTO 100
855
          ENDIF
856
          IF (sigma /= sigma) THEN
857
            CALL Fatal( 'RealBiCGStab(l)', 'Breakdown error: sigma == NaN.' )
858
          ENDIF
859
          
860
          alpha = rho1/sigma
861

862
          !$OMP PARALLEL PRIVATE(j)
863
          !$OMP DO SCHEDULE(STATIC)
864
          DO i=1,n
865
             x(i) = x(i) + alpha * work(i,u)
866
          END DO
867
          !$OMP END DO NOWAIT
868
          DO j=0,k-1
869
             !$OMP DO SCHEDULE(STATIC)
870
             DO i=1,n
871
                work(i,r+j) = work(i,r+j) - alpha * work(i,u+j+1)
872
             END DO
873
             !$OMP END DO
874
          ENDDO
875
          !$OMP END PARALLEL
876

877
          ! CALL C_rpcond( t, work(:,r+k-1), ipar,pcondrsubr )
878
          ! CALL C_matvec( t, work(:,r+k), ipar, matvecsubr )
879
          CALL C_rpcond( t, work(1,r+k-1), ipar,pcondrsubr )
880
          CALL C_matvec( t, work(1,r+k), ipar, matvecsubr )
881

882
          ! rnrm = normfun(n, work(1:n,r), 1 )
883
          rnrm = normfun(n, work(1,r), 1 )
884

885
          IF (rnrm /= rnrm) THEN
886
            CALL Fatal( 'RealBiCGStab(l)', 'Breakdown error: rnrm == NaN.' )
887
          ENDIF
888

889
          mxnrmx = MAX (mxnrmx, rnrm)
890
          mxnrmr = MAX (mxnrmr, rnrm)
891

892
          !----------------------------------------------------------------------
893
          ! In some simple cases, a few BiCG updates may already be enough to
894
          ! obtain the solution. The following is for handling this special case. 
895
          !----------------------------------------------------------------------
896
          errorind = rnrm / bnrm
897
          
898
!         IF( OutputInterval /= 0) THEN
899
!           WRITE (*, '(I8, 2E11.4)') 0, rnrm, errorind
900
!         END IF
901

902
          Converged = (errorind < Tol) 
903
          Diverged = (errorind /= errorind)
904

905
          IF (Converged .OR. Diverged) THEN
906
             EarlyExit = .TRUE.
907
             EXIT
908
          END IF
909
        END DO
910

911
        IF (EarlyExit) EXIT
912

913
        !--------------------------------------
914
        ! --- The convex polynomial part ---
915
        !--------------------------------------
916
        DO i=1,l+1
917
          DO j=1,i
918
             ! rwork(i,j) = dotprodfun(n, work(1:n,r+i-1), 1, work(1:n,r+j-1), 1 ) 
919
             rwork(i,j) = dotprodfun(n, work(1,r+i-1), 1, work(1,r+j-1), 1 ) 
920
          END DO
921
        END DO
922
        DO j=2,l+1
923
           DO i=1,j-1
924
              rwork(i,j) = rwork(j,i)
925
           END DO
926
        END DO
927
        DO j=0,l-1
928
           DO i=1,l+1
929
              rwork(i,zz+j) = rwork(i,z+j)
930
           END DO
931
        END DO
932
        DO j=1,l-1
933
           DO i=1,l-1
934
              tmpmtr(i,j) = rwork(i+1,zz+j)
935
           END DO
936
        END DO
937
        ! CALL dgetrf (l-1, l-1, rwork(2:l,zz+1:zz+l-1), l-1, &
938
        !     iwork, stat)
939
        CALL dgetrf (l-1, l-1, tmpmtr, l-1, &
940
             iwork, stat)
941
      
942
        ! --- tilde r0 and tilde rl (small vectors)
943
        
944
        rwork(1,y0) = -one
945
        DO i=2,l
946
           rwork(i,y0) = rwork(i,z)
947
        END DO
948
        DO i=1,l-1
949
           tmpvec(i) = rwork(i+1,y0)
950
        END DO
951
        ! CALL dgetrs('n', l-1, 1, rwork(2:l,zz+1:zz+l-1), l-1, iwork, &
952
        !     rwork(2:l,y0), l-1, stat)
953
        CALL dgetrs('n', l-1, 1, tmpmtr, l-1, iwork, &
954
             tmpvec, l-1, stat)
955
        DO i=1,l-1
956
           rwork(i+1,y0) = tmpvec(i)
957
        END DO
958
        rwork(l+1,y0) = zero
959
        
960
        rwork(1,yl) = zero
961
        DO i=1,l-1
962
           rwork(i+1,yl) = rwork(i+1,z+l)
963
           tmpvec(i) = rwork(i+1,yl)
964
        END DO
965
        ! CALL dgetrs ('n', l-1, 1, rwork(2:l,zz+1:zz+l-1), l-1, iwork, &
966
        !     rwork(2:l,yl), l-1, stat)
967
        CALL dgetrs ('n', l-1, 1, tmpmtr, l-1, iwork, &
968
             tmpvec, l-1, stat)
969
        DO i=1,l-1
970
           rwork(i+1,yl) = tmpvec(i)
971
        END DO
972
        rwork(l+1,yl) = -one
973
      
974
        ! --- Convex combination          
975
        
976
        CALL dsymv ('u', l+1, one, rwork(1,z), l+1, &
977
            rwork(1,y0), 1, zero, rwork(1,y), 1)
978
        kappa0 = ddot(l+1, rwork(1,y0), 1, rwork(1,y), 1)
979

980
        ! If untreated this would result to NaN's
981
        IF( kappa0 <= 0.0 ) THEN
982
          CALL Warn('RealBiCGStab(l)','kappa0^2 is non-positive, iteration halted')
983
          Halted = .TRUE.
984
          GOTO 100
985
        END IF
986
        kappa0 = SQRT( kappa0 ) 
987

988
        CALL dsymv ('u', l+1, one, rwork(1,z), l+1, &
989
            rwork(1,yl), 1, zero, rwork(1,y), 1)
990
        kappal = ddot(l+1, rwork(1,yl), 1, rwork(1,y), 1 )
991
        
992
        ! If untreated this would result to NaN's
993
        IF( kappal <= 0.0 ) THEN
994
          CALL Warn('RealBiCGStab(l)','kappal^2 is non-positive, iteration halted')
995
          Halted = .TRUE.
996
          GOTO 100 
997
        END IF
998
        kappal = SQRT( kappal )
999

1000
        CALL dsymv ('u', l+1, one, rwork(1,z), l+1, &
1001
          rwork(1,y0), 1, zero, rwork(1,y), 1)
1002

1003
        varrho = ddot(l+1, rwork(1,yl), 1, rwork(1,y), 1) / &
1004
            (kappa0*kappal)
1005
        
1006
        hatgamma = varrho/ABS(varrho) * MAX(ABS(varrho),7d-1) * &
1007
            kappa0/kappal
1008
        DO i=1,l+1
1009
           rwork(i,y0) = rwork(i,y0) - hatgamma * rwork(i,yl)
1010
        END DO
1011
        !  --- Update
1012
        
1013
        omega = rwork(l+1,y0)
1014
        !$OMP PARALLEL PRIVATE(j,i) FIRSTPRIVATE(rwork)
1015
        DO j=1,l
1016
           !$OMP DO SCHEDULE(STATIC)
1017
           DO i=1,n
1018
              work(i,u) = work(i,u) - rwork(j+1,y0) * work(i,u+j)
1019
           END DO
1020
           !$OMP END DO
1021
           !$OMP DO SCHEDULE(STATIC)
1022
           DO i=1,n
1023
              x(i) = x(i) + rwork(j+1,y0) * work(i,r+j-1)
1024
           END DO
1025
           !$OMP END DO
1026
           !$OMP DO SCHEDULE(STATIC)
1027
           DO i=1,n
1028
              work(i,r) = work(i,r) - rwork(j+1,y0) * work(i,r+j)
1029
           END DO
1030
           !$OMP END DO
1031
        ENDDO
1032
        !$OMP END PARALLEL
1033
    
1034
        CALL dsymv ('u', l+1, one, rwork(1,z), l+1, &
1035
            rwork(1,y0), 1, zero, rwork(1,y), 1)
1036
        rnrm = ddot(l+1, rwork(1,y0), 1, rwork(1,y), 1)
1037

1038
        IF( rnrm < 0.0 ) THEN
1039
          CALL Warn('RealBiCGStab(l)','rnrm^2 is negative, iteration halted')
1040
          Halted = .TRUE.
1041
          GOTO 100 
1042
        END IF        
1043
        rnrm = SQRT( rnrm ) 
1044
        
1045
        !---------------------------------------
1046
        !  --- The reliable update part ---
1047
        !---------------------------------------
1048
        
1049
        mxnrmx = MAX (mxnrmx, rnrm)
1050
        mxnrmr = MAX (mxnrmr, rnrm)
1051
        xpdt = (rnrm < delta*rnrm0 .AND. rnrm0 < mxnrmx)
1052
        rcmp = ((rnrm < delta*mxnrmr .AND. rnrm0 < mxnrmr) .OR. xpdt)
1053
        IF (rcmp) THEN
1054
          ! PRINT *, 'Performing residual update...'
1055
          CALL C_rpcond( t, x, ipar,pcondrsubr )
1056
          ! CALL C_matvec( t, work(:,r), ipar, matvecsubr )
1057
          CALL C_matvec( t, work(1,r), ipar, matvecsubr )
1058

1059
          mxnrmr = rnrm
1060
          !$OMP PARALLEL DO SCHEDULE(STATIC)
1061
          DO i=1,n
1062
             work(i,r) = work(i,bp) - work(i,r)
1063
          END DO
1064
          !$OMP END PARALLEL DO
1065
          IF (xpdt) THEN
1066
            ! PRINT *, 'Performing solution update...'
1067
            !$OMP PARALLEL DO SCHEDULE(STATIC)
1068
             DO i=1,n
1069
                work(i,xp) = work(i,xp) + t(i)
1070
                x(i) = zero
1071
                work(i,bp) = work(i,r)
1072
             END DO
1073
             !$OMP END PARALLEL DO
1074

1075
             mxnrmx = rnrm
1076
          ENDIF
1077
        ENDIF
1078
        
1079
        IF (rcmp) THEN
1080
          IF (xpdt) THEN       
1081
             !$OMP PARALLEL DO SCHEDULE(STATIC)
1082
             DO i=1,n
1083
                t(i) = work(i,xp)
1084
             END DO
1085
             !$OMP END PARALLEL DO
1086
          ELSE
1087
             !$OMP PARALLEL DO SCHEDULE(STATIC)
1088
             DO i=1,n
1089
                t(i) = t(i) + work(i,xp)  
1090
             END DO
1091
             !$OMP END PARALLEL DO
1092
          END IF
1093
        ELSE
1094
          CALL C_rpcond( t, x, ipar,pcondrsubr )
1095
          !$OMP PARALLEL DO
1096
          DO i=1,n
1097
             t(i) = t(i)+work(i,xp)
1098
          END DO
1099
          !$OMP END PARALLEL DO
1100
        END IF
1101
      
1102
        errorind = rnrm / bnrm
1103

1104
        IF( MOD(Round,OutputInterval) == 0) THEN
1105
          WRITE (*, '(I8, 2E11.4)') Round, rnrm, errorind
1106
          CALL FLUSH(6)
1107
        END IF
1108
        
1109
        IF( Robust ) THEN
1110
          IF (Round>=RobustStart ) THEN
1111
            IF( errorInd < RobustStep * BestNorm ) THEN
1112
              BestIter = Round
1113
              BestNorm = errorInd
1114
              Bestx = x
1115
              BadIterCount = 0
1116
            ELSE
1117
              BadIterCount = BadIterCount + 1
1118
            END IF
1119

1120
            IF( BestNorm <  RobustTol .AND. &
1121
                  ( errorInd > RobustMaxTol .OR. BadIterCount > MaxBadIter ) ) THEN
1122
              EXIT
1123
            END IF
1124
          END IF
1125
        END IF
1126
               
1127
        Converged = (errorind < Tol) 
1128
        Diverged = (errorind > MaxTol) .OR. (errorind /= errorind)
1129
        IF( Converged .OR. Diverged) EXIT    
1130
      END DO
1131

1132
100   IF( Robust ) THEN
1133
        IF( BestNorm < RobustTol ) THEN
1134
          Converged = .TRUE.
1135
        END IF
1136
        IF( BestNorm < errorInd ) THEN
1137
          x = Bestx
1138
        END IF
1139
        IF(OutputInterval /= HUGE(OutputInterval)) THEN
1140
          WRITE(*,'(A,I8,E11.4,I8,2E11.4)') 'BiCGStabl robust: ',&
1141
              MIN(MaxRounds,Round), BestNorm, BestIter, rnrm, errorind
1142
          CALL FLUSH(6)
1143
        END IF
1144
      ELSE
1145
        IF(OutputInterval /= HUGE(OutputInterval)) THEN
1146
          WRITE (*, '(A, I8, 2E11.4)') 'BiCGStabl: ', MIN(MaxRounds,Round), rnrm, errorind
1147
          CALL FLUSH(6)
1148
        END IF
1149
      END IF
1150
            
1151
      !------------------------------------------------------------
1152
      ! We have solved z = P*x, with P the preconditioner, so finally 
1153
      ! solve the true unknown x
1154
      !------------------------------------------------------------
1155
      !$OMP PARALLEL DO
1156
      DO i=1,n
1157
         t(i) = x(i)
1158
      END DO
1159
      !$OMP END PARALLEL DO
1160
      CALL C_rpcond( x, t, ipar,pcondrsubr )
1161
      !$OMP PARALLEL DO
1162
      DO i=1,n
1163
         x(i) = x(i) + work(i,xp)
1164
      END DO
1165
      !$OMP END PARALLEL DO
1166
!------------------------------------------------------------------------------
1167
    END SUBROUTINE RealBiCGStabl
1168
!------------------------------------------------------------------------------
1169

1170
  END SUBROUTINE itermethod_bicgstabl
1171
!------------------------------------------------------------------------------
1172

1173

1174

1175
!------------------------------------------------------------------------------
1176
!>   This routine solves real linear systems Ax = b by using the GCR algorithm 
1177
!> (Generalized Conjugate Residual).
1178
!------------------------------------------------------------------------------
1179
 SUBROUTINE itermethod_gcr( xvec, rhsvec, &
1180
      ipar, dpar, work, matvecsubr, pcondlsubr, &
1181
      pcondrsubr, dotprodfun, normfun, stopcfun )
1182

1183
    USE huti_interfaces
1184
    IMPLICIT NONE
1185
    PROCEDURE( mv_iface_d ), POINTER :: matvecsubr
1186
    PROCEDURE( pc_iface_d ), POINTER :: pcondlsubr
1187
    PROCEDURE( pc_iface_d ), POINTER :: pcondrsubr
1188
    PROCEDURE( dotp_iface_d ), POINTER :: dotprodfun
1189
    PROCEDURE( norm_iface_d ), POINTER :: normfun
1190
    PROCEDURE( stopc_iface_d ), POINTER :: stopcfun
1191

1192
    INTEGER, DIMENSION(HUTI_IPAR_DFLTSIZE) :: ipar
1193
    REAL(KIND=dp), TARGET, DIMENSION(HUTI_NDIM) :: xvec, rhsvec
1194
    DOUBLE PRECISION, DIMENSION(HUTI_DPAR_DFLTSIZE) :: dpar
1195
    REAL(KIND=dp), DIMENSION(HUTI_WRKDIM,HUTI_NDIM) :: work
1196

1197
    INTEGER :: ndim, RestartN
1198
    INTEGER :: Rounds, MinIter, OutputInterval
1199
    REAL(KIND=dp) :: MinTol, MaxTol, Residual
1200
    LOGICAL :: Converged, Diverged, UseStopCFun
1201
    LOGICAL :: PseudoComplex
1202

1203
    TYPE(Matrix_t),POINTER::A
1204

1205
    REAL(KIND=dp), POINTER :: x(:),b(:)
1206

1207
    CALL Info('Itermethod_gcr','Starting GCR iteration',Level=25)
1208

1209
    
1210
    ndim = HUTI_NDIM
1211
    Rounds = HUTI_MAXIT
1212
    MinIter = HUTI_MINIT
1213
    MinTol = HUTI_TOLERANCE
1214
    MaxTol = HUTI_MAXTOLERANCE
1215
    OutputInterval = HUTI_DBUGLVL
1216
    RestartN = HUTI_GCR_RESTART 
1217
    UseStopCFun = HUTI_STOPC == HUTI_USUPPLIED_STOPC
1218

1219
    Converged = .FALSE.
1220
    Diverged = .FALSE.
1221
    PseudoComplex = ( HUTI_PSEUDOCOMPLEX > 0 )  
1222
      
1223
    x => xvec
1224
    b => rhsvec
1225
    nc = 0
1226

1227
    A => GlobalMatrix
1228
    CM => A % ConstraintMatrix
1229
    Constrained = ASSOCIATED(CM)
1230
    
1231
    IF (Constrained) THEN
1232
      nc = CM % NumberOfRows
1233
      Constrained = nc>0
1234
      IF(Constrained) THEN
1235
        ALLOCATE(x(ndim+nc),b(ndim+nc))
1236
        IF(.NOT.ALLOCATED(CM % ExtraVals))THEN
1237
          ALLOCATE(CM % ExtraVals(nc)); CM % extraVals=0._dp
1238
        END IF
1239
        b(1:ndim) = rhsvec; b(ndim+1:) = CM % RHS
1240
        x(1:ndim) = xvec; x(ndim+1:) = CM % extraVals
1241
      END IF
1242
    END IF
1243
    
1244
    CALL GCR(ndim+nc, GlobalMatrix, x, b, Rounds, MinTol, MaxTol, Residual, &
1245
        Converged, Diverged, OutputInterval, RestartN, MinIter )
1246

1247
    
1248
    IF(Constrained) THEN
1249
      xvec = x(1:ndim)
1250
      rhsvec = b(1:ndim)
1251
      CM % extraVals = x(ndim+1:ndim+nc)
1252
      DEALLOCATE(x,b)
1253
    END IF
1254

1255
    IF(Converged) HUTI_INFO = HUTI_CONVERGENCE
1256
    IF(Diverged) HUTI_INFO = HUTI_DIVERGENCE
1257
    IF ( (.NOT. Converged) .AND. (.NOT. Diverged) ) HUTI_INFO = HUTI_MAXITER   
1258

1259
  CONTAINS 
1260
    
1261
    
1262
    SUBROUTINE GCR( n, A, x, b, Rounds, MinTolerance, MaxTolerance, Residual, &
1263
        Converged, Diverged, OutputInterval, m, MinIter) 
1264
!------------------------------------------------------------------------------
1265
      TYPE(Matrix_t), POINTER :: A
1266
      INTEGER :: Rounds, MinIter
1267
      REAL(KIND=dp) :: x(n),b(n)
1268
      LOGICAL :: Converged, Diverged
1269
      REAL(KIND=dp) :: MinTolerance, MaxTolerance, Residual
1270
      INTEGER :: n, OutputInterval, m
1271
      REAL(KIND=dp) :: bnorm,rnorm
1272
      REAL(KIND=dp), ALLOCATABLE :: R(:)
1273

1274
      REAL(KIND=dp), ALLOCATABLE :: S(:,:), V(:,:), T1(:), T2(:)
1275

1276
!------------------------------------------------------------------------------
1277
      INTEGER :: i,j,k,ksum=0
1278
      REAL(KIND=dp) :: alpha, beta, trueres(n), trueresnorm, normerr
1279
      REAL(KIND=dp) :: beta_im
1280
!------------------------------------------------------------------------------
1281
      INTEGER :: allocstat
1282
        
1283
      ALLOCATE( R(n), T1(n), T2(n), STAT=allocstat )
1284
      IF( allocstat /= 0 ) THEN
1285
        CALL Fatal('GCR','Failed to allocate memory of size: '//I2S(n))
1286
      END IF
1287

1288
      IF ( m > 1 ) THEN
1289
        ALLOCATE( S(n,m-1), V(n,m-1), STAT=allocstat )
1290
        IF( allocstat /= 0 ) THEN
1291
          CALL Fatal('GCR','Failed to allocate memory of size: '&
1292
              //I2S(n)//' x '//I2S(m-1))
1293
        END IF
1294
        
1295
         V(1:n,1:m-1) = 0.0d0	
1296
         S(1:n,1:m-1) = 0.0d0
1297
      END IF	
1298
      
1299
      CALL C_matvec( x, r, ipar, matvecsubr )
1300
      r(1:n) = b(1:n) - r(1:n)
1301
      
1302
      bnorm = normfun(n, b, 1)
1303
      rnorm = normfun(n, r, 1)
1304
      
1305
      IF (UseStopCFun) THEN
1306
        Residual = stopcfun(x,b,r,ipar,dpar)
1307
      ELSE
1308
        Residual = rnorm / bnorm
1309
      END IF
1310
      Converged = (Residual < MinTolerance) .AND. ( MinIter <= 0 )
1311
      Diverged = (Residual > MaxTolerance) .OR. (Residual /= Residual)
1312
      IF( Converged .OR. Diverged) RETURN
1313
      
1314
      DO k=1,Rounds
1315
         !----------------------------------------------
1316
         ! Check for restarting
1317
         !----------------------------------------------
1318
         IF ( MOD(k,m)==0 ) THEN
1319
            j = m
1320
         ELSE
1321
            j = MOD(k,m)
1322
            !--------------------------------------------
1323
            ! Compute the true residual when restarting:
1324
            !--------------------------------------------
1325
            IF ( (j==1) .AND. (k>1) ) THEN
1326
               CALL C_matvec( x, r, ipar, matvecsubr )
1327
               r(1:n) = b(1:n) - r(1:n)
1328
            END IF
1329
         END IF
1330

1331
         !----------------------------------------------------------
1332
         ! Perform the preconditioning...
1333
         !---------------------------------------------------------------
1334
         CALL C_rpcond( T1, r, ipar, pcondrsubr )
1335
         CALL C_matvec( T1, T2, ipar, matvecsubr )
1336

1337
         !--------------------------------------------------------------
1338
         ! Perform the orthogonalization of the search directions....
1339
         !--------------------------------------------------------------
1340
         DO i=1,j-1
1341
           beta = dotprodfun(n, V(1:n,i), 1, T2(1:n), 1 )
1342

1343
           IF( PseudoComplex ) THEN
1344
             ! The even call is for the complex part of beta
1345
             ! This has to be before the T1 and T2 vectors are tampered
1346
             ! For convenience we subtract 
1347
             beta_im = dotprodfun(n, V(1:n,i), 1, T2(1:n), 1 )
1348

1349
             IF( HUTI_PSEUDOCOMPLEX == 2 ) THEN
1350
               T1(1:n/2) = T1(1:n/2) + beta_im * S(n/2+1:n,i) 
1351
               T1(n/2+1:n) = T1(n/2+1:n) - beta_im * S(1:n/2,i)                    
1352
               
1353
               T2(1:n/2) = T2(1:n/2) + beta_im * V(1+n/2:n,i)
1354
               T2(1+n/2:n) = T2(1+n/2:n) - beta_im * V(1:n/2,i)                                
1355
             ELSE
1356
               T1(1:n:2) = T1(1:n:2) + beta_im * S(2:n:2,i) 
1357
               T1(2:n:2) = T1(2:n:2) - beta_im * S(1:n:2,i)                    
1358
               
1359
               T2(1:n:2) = T2(1:n:2) + beta_im * V(2:n:2,i)
1360
               T2(2:n:2) = T2(2:n:2) - beta_im * V(1:n:2,i)
1361
             END IF
1362
           END IF
1363
           
1364
           T1(1:n) = T1(1:n) - beta * S(1:n,i)
1365
           T2(1:n) = T2(1:n) - beta * V(1:n,i)        
1366
         END DO
1367

1368
         alpha = normfun(n, T2(1:n), 1 )
1369
         T1(1:n) = 1.0d0/alpha * T1(1:n)
1370
         T2(1:n) = 1.0d0/alpha * T2(1:n)
1371

1372
         !-------------------------------------------------------------
1373
         ! The update of the solution and save the search data...
1374
         !------------------------------------------------------------- 
1375
         beta = dotprodfun(n, T2(1:n), 1, r(1:n), 1 )
1376

1377
         IF( PseudoComplex ) THEN
1378
           beta_im = dotprodfun(n, T2(1:n), 1, r(1:n), 1 )
1379

1380
           IF( HUTI_PSEUDOCOMPLEX == 2 ) THEN
1381
             x(1:n/2) = x(1:n/2) - beta_im * T1(1+n/2:n)
1382
             x(1+n/2:n) = x(1+n/2:n) + beta_im * T1(1:n/2)                    
1383
             r(1:n/2) = r(1:n/2) + beta_im * T2(1+n/2:n)
1384
             r(1+n/2:n) = r(1+n/2:n) - beta_im * T2(1:n/2)                    
1385
           ELSE
1386
             x(1:n:2) = x(1:n:2) - beta_im * T1(2:n:2)
1387
             x(2:n:2) = x(2:n:2) + beta_im * T1(1:n:2)                    
1388
             r(1:n:2) = r(1:n:2) + beta_im * T2(2:n:2)
1389
             r(2:n:2) = r(2:n:2) - beta_im * T2(1:n:2)                    
1390
           END IF
1391
         END IF
1392
                      
1393
         x(1:n) = x(1:n) + beta * T1(1:n)      
1394
         r(1:n) = r(1:n) - beta * T2(1:n)
1395

1396
         IF ( j /= m ) THEN
1397
           S(1:n,j) = T1(1:n)
1398
           V(1:n,j) = T2(1:n)
1399
	 END IF       
1400

1401
         !--------------------------------------------------------------
1402
         ! Check whether the convergence criterion is met 
1403
         !--------------------------------------------------------------
1404
         rnorm = normfun(n, r, 1)
1405

1406
         IF (UseStopCFun) THEN
1407
           Residual = stopcfun(x,b,r,ipar,dpar)
1408
           IF( MOD(k,OutputInterval) == 0) THEN
1409
             WRITE (*, '(A, I6, 2E12.4)') '   gcr:',k, rnorm / bnorm, residual
1410
             CALL FLUSH(6)
1411
           END IF
1412
         ELSE
1413
           Residual = rnorm / bnorm
1414
           IF( MOD(k,OutputInterval) == 0) THEN
1415
             IF( PseudoComplex ) THEN
1416
               WRITE (*, '(A, I6, 3E12.4, A)') '   gcr:',k, residual, beta, beta_im,'i'
1417
             ELSE
1418
               WRITE (*, '(A, I6, 2E12.4)') '   gcr:',k, residual, beta
1419
             END IF
1420
             CALL FLUSH(6)
1421
           END IF
1422
         END IF
1423
           
1424
         Converged = (Residual < MinTolerance) .AND. ( k >= MinIter )
1425
         !-----------------------------------------------------------------
1426
         ! Make an additional check that the true residual agrees with 
1427
         ! the iterated residual:
1428
         !-----------------------------------------------------------------
1429
         IF (Converged ) THEN
1430
           CALL C_matvec( x, trueres, ipar, matvecsubr )
1431
           trueres(1:n) = b(1:n) - trueres(1:n)
1432
           TrueResNorm = normfun(n, trueres, 1)
1433
           NormErr = ABS(TrueResNorm - rnorm)/TrueResNorm
1434
           
1435
           IF ( NormErr > 1.0d-1 ) THEN
1436
             CALL Warn('IterMethod_GCR','Iterated GCR solution may not be accurate')
1437
             i = 4
1438
           ELSE
1439
             i = 8
1440
           END IF
1441
           WRITE( Message,'(A,I0,A,ES12.3)') 'Iterated residual norm after ',k,' iters:', rnorm
1442
           CALL Info('IterMethod_GCR', Message, Level=i)
1443
           WRITE( Message,'(A,ES12.3)') 'True residual norm::', TrueResNOrm
1444
           CALL Info('IterMethod_GCR', Message, Level=i)            
1445

1446
           IF( InfoActive(20) ) THEN
1447
             ksum = ksum + k
1448
             CALL Info('IterMethod_GCR','Total number of GCR iterations: '//I2S(ksum))           
1449
           END IF
1450
           
1451
         END IF
1452
         Diverged = (Residual > MaxTolerance) .OR. (Residual /= Residual)    
1453
         IF( Converged .OR. Diverged) EXIT
1454
        
1455
      END DO
1456
      
1457
      DEALLOCATE( R, T1, T2 )
1458
      IF ( m > 1 ) DEALLOCATE( S, V)
1459
      
1460
    END SUBROUTINE GCR
1461
    
1462
!------------------------------------------------------------------------------
1463
  END SUBROUTINE itermethod_gcr
1464
!------------------------------------------------------------------------------
1465

1466

1467
   
1468
!-----------------------------------------------------------------------------------
1469
!>  This subroutine solves real linear systems Ax = b by using the IDR(s) algorithm
1470
!>  with s >= 1 and the right-oriented preconditioning.
1471
!------------------------------------------------------------------------------
1472
  SUBROUTINE itermethod_idrs( xvec, rhsvec, &
1473
      ipar, dpar, work, matvecsubr, pcondlsubr, &
1474
      pcondrsubr, dotprodfun, normfun, stopcfun )
1475
!------------------------------------------------------------------------------
1476
    USE huti_interfaces
1477
    IMPLICIT NONE
1478
    PROCEDURE( mv_iface_d ), POINTER :: matvecsubr
1479
    PROCEDURE( pc_iface_d ), POINTER :: pcondlsubr
1480
    PROCEDURE( pc_iface_d ), POINTER :: pcondrsubr
1481
    PROCEDURE( dotp_iface_d ), POINTER :: dotprodfun
1482
    PROCEDURE( norm_iface_d ), POINTER :: normfun
1483
    PROCEDURE( stopc_iface_d ), POINTER :: stopcfun
1484

1485
    INTEGER, DIMENSION(HUTI_IPAR_DFLTSIZE) :: ipar
1486
    REAL(KIND=dp), TARGET, DIMENSION(HUTI_NDIM) :: xvec, rhsvec
1487
    ! DOUBLE PRECISION, DIMENSION(HUTI_NDIM), TARGET :: xvec, rhsvec
1488
    DOUBLE PRECISION, DIMENSION(HUTI_DPAR_DFLTSIZE) :: dpar
1489
    REAL(KIND=dp), DIMENSION(HUTI_WRKDIM,HUTI_NDIM) :: work
1490
    ! DOUBLE PRECISION, DIMENSION(HUTI_WRKDIM,HUTI_NDIM) :: work
1491
    INTEGER :: ndim,i,j,k
1492
    INTEGER :: Rounds, OutputInterval, s
1493
    REAL(KIND=dp) :: MinTol, MaxTol, Residual
1494
    LOGICAL :: Converged, Diverged, UseStopCFun
1495

1496
    TYPE(Matrix_t), POINTER :: A
1497

1498
    REAL(KIND=dp), POINTER :: x(:),b(:)
1499

1500
    ! Variables related to robust mode
1501
    LOGICAL :: Robust 
1502
    INTEGER :: BestIter,BadIterCount,MaxBadIter
1503
    REAL(KIND=dp) :: BestNorm,RobustStep,RobustTol,RobustMaxTol
1504
    REAL(KIND=dp), ALLOCATABLE :: Bestx(:)
1505

1506
    LOGICAL :: Smoothing 
1507

1508
    A => GlobalMatrix
1509
    CM => A % ConstraintMatrix
1510
    Constrained = ASSOCIATED(CM)
1511
    
1512
    ndim = HUTI_NDIM
1513

1514
    x => xvec
1515
    b => rhsvec
1516
    nc = 0
1517
    IF (Constrained) THEN
1518
      nc = CM % NumberOfRows
1519
      Constrained = nc>0
1520
      IF(Constrained) THEN
1521
        ALLOCATE(x(ndim+nc),b(ndim+nc))
1522
        IF(.NOT.ALLOCATED(CM % ExtraVals))THEN
1523
          ALLOCATE(CM % ExtraVals(nc)); CM % extraVals=0._dp
1524
        END IF
1525
        b(1:ndim) = rhsvec; b(ndim+1:) = CM % RHS
1526
        x(1:ndim) = xvec; x(ndim+1:) = CM % extraVals
1527
      END IF
1528
    END IF
1529

1530
    Rounds = HUTI_MAXIT
1531
    MinTol = HUTI_TOLERANCE
1532
    MaxTol = HUTI_MAXTOLERANCE
1533
    OutputInterval = HUTI_DBUGLVL 
1534
    s = HUTI_IDRS_S
1535
    UseStopCFun = HUTI_STOPC == HUTI_USUPPLIED_STOPC
1536
    
1537
    Robust = ( HUTI_ROBUST == 1 )
1538
    IF( Robust ) THEN
1539
      RobustTol = HUTI_ROBUST_TOLERANCE
1540
      RobustStep = HUTI_ROBUST_STEPSIZE
1541
      RobustMaxTol = HUTI_ROBUST_MAXTOLERANCE
1542
      MaxBadIter = HUTI_ROBUST_MAXBADIT
1543
      BestNorm = SQRT(HUGE(BestNorm))
1544
      BadIterCount = 0
1545
      BestIter = 0      
1546
      ALLOCATE( BestX(ndim))
1547
    END IF
1548

1549
    Smoothing = ( HUTI_SMOOTHING == 1) 
1550

1551
    Converged = .FALSE.
1552
    Diverged = .FALSE.
1553
    
1554
    CALL RealIDRS(ndim+nc, A,x,b, Rounds, MinTol, MaxTol, &
1555
         Converged, Diverged, OutputInterval, s )
1556

1557
    IF(Constrained) THEN
1558
      xvec=x(1:ndim)
1559
      rhsvec=b(1:ndim)
1560
      CM % extraVals = x(ndim+1:ndim+nc)
1561
      DEALLOCATE(x,b)
1562
    END IF
1563

1564
    IF( Robust ) THEN
1565
      DEALLOCATE( BestX )
1566
    END IF
1567
    
1568

1569
    IF(Converged) HUTI_INFO = HUTI_CONVERGENCE
1570
    IF(Diverged) HUTI_INFO = HUTI_DIVERGENCE
1571
    IF ( (.NOT. Converged) .AND. (.NOT. Diverged) ) HUTI_INFO = HUTI_MAXITER
1572

1573
  CONTAINS
1574

1575
!-----------------------------------------------------------------------------------
1576
!   The subroutine RealIDRS solves real linear systems Ax = b by using the IDR(s) 
1577
!   algorithm with s >= 1 and the right-oriented preconditioning.
1578
!
1579
!   The subroutine RealIDRS has been written by M.B. van Gijzen
1580
!----------------------------------------------------------------------------------- 
1581
    SUBROUTINE RealIDRS( n, A, x, b, MaxRounds, Tol, MaxTol, Converged, &
1582
        Diverged, OutputInterval, s)
1583
!----------------------------------------------------------------------------------- 
1584
      INTEGER :: s   ! IDR parameter
1585
      INTEGER :: n, MaxRounds, OutputInterval   
1586
      LOGICAL :: Converged, Diverged
1587
      TYPE(Matrix_t), POINTER :: A
1588
      REAL(KIND=dp) :: x(n), b(n)
1589
      REAL(KIND=dp) :: Tol, MaxTol
1590

1591
      ! Local arrays:
1592
      REAL(kind=dp) :: P(n,s)
1593
      REAL(kind=dp) :: G(n,s)
1594
      REAL(kind=dp) :: U(n,s)
1595
      REAL(kind=dp) :: r(n)
1596
      REAL(kind=dp) :: v(n)
1597
      REAL(kind=dp) :: t(n)
1598
      REAL(kind=dp) :: M(s,s), f(s), mu(s)
1599
      REAL(kind=dp) :: alpha(s), beta(s), gamma(s)
1600

1601
      REAL(kind=dp) :: om, tr, tr_s, tt
1602
      REAL(kind=dp) :: nr, nt, rho, kappa
1603

1604
      REAL(kind=dp), ALLOCATABLE :: r_s(:), x_s(:)
1605
      REAL(kind=dp) :: theta
1606
      
1607
      INTEGER :: iter                         ! number of iterations
1608
      INTEGER :: ii                           ! inner iterations index
1609
      INTEGER :: jj                           ! G-space index
1610
      REAL(kind=dp) :: normb, normr, errorind ! for tolerance check
1611
      INTEGER :: i,j,k,l                      ! loop counters
1612
!----------------------------------------------------------------------------------- 
1613
      
1614
      ! Compute initial residual
1615
      normb = normfun(n,b,1)
1616
      CALL C_matvec( x, t, ipar, matvecsubr )
1617
      r = b - t
1618
      
1619
      !-------------------------------------------------------------------
1620
      ! Check whether the initial guess satisfies the stopping criterion
1621
      !--------------------------------------------------------------------
1622
      IF (UseStopCFun) THEN
1623
        errorind = stopcfun(x,b,r,ipar,dpar)
1624
      ELSE
1625
        normr = normfun(n,r,1)
1626
        errorind = normr / normb
1627
      END IF
1628
      Converged = (errorind < Tol)
1629
      Diverged = (errorind > MaxTol) .OR. (errorind /= errorind)
1630

1631
      IF ( Converged .OR. Diverged ) RETURN
1632

1633
      U = 0.0d0
1634
      IF ( Smoothing ) THEN
1635
        ALLOCATE( r_s(n), x_s(n) )
1636
        x_s = x
1637
        r_s = r 
1638
      END IF
1639

1640
      ! Define P(n,s) and kappa
1641
#if 1
1642
      CALL RANDOM_NUMBER(P)
1643
#else
1644
      ! this is alternative generation of initial basis vectors
1645
      ! it is deterministic but not as good...
1646
      l = 0
1647
      k = 2        
1648
      DO j=1,s
1649
        DO i=1,n
1650
          P(i,j) = MODULO(i+l,k) / (1.0*(k-1)) 
1651
        END DO
1652
        l = k
1653
        k = 2*k + 1
1654
      END DO
1655
#endif
1656
              
1657
      DO j = 1,s
1658
        DO k = 1,j-1
1659
          alpha(k) = dotprodfun(n, P(:,k), 1, P(:,j), 1 )
1660
          P(:,j) = P(:,j) - alpha(k)*P(:,k)
1661
        END DO
1662
        P(:,j) = P(:,j)/normfun(n,P(:,j),1)
1663
      END DO
1664
      kappa = 0.7d0
1665

1666
      ! Initialize local variables:
1667
      M = 0.0d0
1668
      om = 1.0d0
1669
      iter = 0
1670
      jj = 0
1671
      ii = 0
1672

1673
      
1674
      ! This concludes the initialisation phase    
1675
      
1676
      ! Main iteration loop, build G-spaces:
1677
      DO WHILE ( (.NOT. Converged) .AND. (.NOT. Diverged) ) 
1678

1679
        !!+++++++++++++++++++++++++++++++++++++++++++++++++++++++
1680
        ! Generate s vectors in G_j
1681
        !!+++++++++++++++++++++++++++++++++++++++++++++++++++++++
1682

1683
        ! New right-hand side for small system:
1684
        DO k = 1,s
1685
          f(k) = dotprodfun(n, P(:,k), 1, r, 1 )
1686
        END DO
1687

1688
        DO k = 1,s
1689

1690
          ! Update inner iteration counter
1691
          ii = ii + 1
1692

1693
          ! Compute new v
1694
          v = r
1695
          IF ( jj > 0 ) THEN
1696

1697
            ! Solve small system (Note: M is lower triangular) and make v orthogonal to P:
1698
            DO i = k,s
1699
              gamma(i) = f(i)
1700
              DO j = k,i-1
1701
                gamma(i) = gamma(i) - M(i,j)*gamma(j)
1702
              END DO
1703
              gamma(i) = gamma(i)/M(i,i)
1704
              v = v - gamma(i)*G(:,i)
1705
            END DO
1706

1707
            ! Compute new U(:,k)
1708
            CALL C_rpcond( t, v, ipar, pcondrsubr ) 
1709
            t = om*t
1710
            DO i = k,s
1711
              t = t + gamma(i)*U(:,i)
1712
            END DO
1713
            U(:,k) = t
1714

1715
          ELSE
1716

1717
            ! Updates for the first s iterations (in G_0):
1718
            CALL C_rpcond( U(:,k), v, ipar, pcondrsubr )
1719

1720
          END IF
1721

1722
          ! Compute new G(:,k), G(:,k) is in space G_j
1723
          CALL C_matvec( U(:,k), G(:,k), ipar, matvecsubr )
1724

1725
          ! Bi-Orthogonalise the new basis vectors:
1726
          DO i = 1,s
1727
            mu(i) = dotprodfun(n, P(:,i), 1, G(:,k), 1 )
1728
          END DO
1729
          DO i = 1,k-1
1730
            alpha(i) = mu(i)
1731
            DO j = 1, i-1
1732
              alpha(i) = alpha(i) - M(i,j)*alpha(j)
1733
            END DO
1734
            alpha(i) = alpha(i)/M(i,i)
1735
            G(:,k) = G(:,k) - G(:,i)*alpha(i)
1736
            U(:,k) = U(:,k) - U(:,i)*alpha(i)
1737
            mu(k:s)  = mu(k:s)  - M(k:s,i)*alpha(i)
1738
          END DO
1739
          M(k:s,k) = mu(k:s)
1740

1741
          ! Break down?
1742
          IF ( ABS(M(k,k)) <= TINY(tol) ) THEN
1743
            Diverged = .TRUE.
1744
            EXIT
1745
          END IF
1746

1747
          ! Make r orthogonal to p_i, i = 1..k, update solution and residual
1748
          beta(k) = f(k)/M(k,k)
1749
          r = r - beta(k)*G(:,k)
1750
          x = x + beta(k)*U(:,k)
1751
          
1752
          ! New f = P'*r (first k  components are zero)
1753
          IF ( k < s ) THEN
1754
            f(k+1:s)   = f(k+1:s) - beta(k)*M(k+1:s,k)
1755
          END IF
1756

1757
          IF (Smoothing) THEN
1758
            t = r_s - r
1759
            tr_s = dotprodfun(n, t, 1, r_s, 1 )
1760
            tt = dotprodfun(n, t, 1, t, 1 )
1761
            theta = tr_s / tt
1762
            
1763
            r_s = r_s - theta * t
1764
            x_s = x_s - theta * (x_s - x)
1765
          END IF
1766

1767
          ! Check for convergence
1768
          iter = iter + 1
1769
          IF (UseStopCFun) THEN
1770
            IF (Smoothing) THEN
1771
              errorind = stopcfun(x_s,b,r_s,ipar,dpar)
1772
            ELSE
1773
              errorind = stopcfun(x,b,r,ipar,dpar)
1774
            END IF
1775
          ELSE
1776
            IF (Smoothing) THEN
1777
              normr = normfun(n,r_s,1)  
1778
            ELSE
1779
              normr = normfun(n,r,1)
1780
            END IF
1781
            errorind = normr/normb
1782
          END IF
1783

1784
          IF( MOD(iter,OutputInterval) == 0) THEN
1785
            WRITE (*, '(I8, E11.4)') iter, errorind
1786
          END IF
1787

1788
          Converged = (errorind < Tol)
1789
          Diverged = (errorind > MaxTol) .OR. (errorind /= errorind)
1790
          IF ( Converged .OR. Diverged ) EXIT
1791
          IF (iter == MaxRounds) EXIT
1792
          
1793
         
1794
        END DO ! Now we have computed s+1 vectors in G_j
1795
        IF ( Converged .OR. Diverged ) EXIT
1796
        IF (iter == MaxRounds) EXIT
1797

1798
        !!+++++++++++++++++++++++++++++++++++++++++++++++++++++++
1799
        ! Compute first residual in G_j+1
1800
        !!+++++++++++++++++++++++++++++++++++++++++++++++++++++++
1801

1802
        ! Update G-space counter
1803
        jj = jj + 1
1804

1805
        ! Compute first residual in G_j+1
1806
        ! Note: r is already perpendicular to P so v = r
1807

1808
        ! Preconditioning:
1809
        CALL C_rpcond( v, r, ipar, pcondrsubr )
1810
        ! Matrix-vector multiplication:
1811
        CALL C_matvec( v, t, ipar, matvecsubr )
1812

1813
        ! Computation of a new omega
1814
        ! 'Maintaining the convergence':
1815
        nr = normfun(n,r,1)
1816
        nt = normfun(n,t,1)
1817
        tr = dotprodfun(n, t, 1, r, 1 )
1818
        rho = ABS(tr/(nt*nr))
1819
        om=tr/(nt*nt)
1820
        IF ( rho < kappa ) THEN
1821
          om = om*kappa/rho
1822
        END IF
1823

1824
        IF ( ABS(om) <= EPSILON(tol) ) THEN
1825
          Diverged = .TRUE.
1826
          EXIT
1827
        END IF
1828

1829
        ! Update solution and residual
1830
        r = r - om*t
1831
        x = x + om*v
1832
        
1833
        IF (Smoothing) THEN
1834
          t = r_s - r
1835
          tr_s = dotprodfun(n, t, 1, r_s, 1 )
1836
          tt = dotprodfun(n, t, 1, t, 1 )
1837
          theta = tr_s / tt
1838
          r_s = r_s - theta * t
1839
          x_s = x_s - theta * (x_s - x)
1840
        END IF
1841

1842
        ! Check for convergence
1843
        iter = iter + 1
1844
        IF (UseStopCFun) THEN
1845
          IF (Smoothing) THEN
1846
            errorind = stopcfun(x_s,b,r_s,ipar,dpar)
1847
          ELSE
1848
            errorind = stopcfun(x,b,r,ipar,dpar)
1849
          END IF
1850
        ELSE
1851
          IF (Smoothing) THEN
1852
            normr = normfun(n,r_s,1)  
1853
          ELSE
1854
            normr = normfun(n,r,1)
1855
          END IF
1856
          errorind = normr/normb
1857
        END IF
1858
        
1859
        IF( MOD(iter,OutputInterval) == 0) THEN
1860
          WRITE (*, '(I8, E11.4)') iter, errorind
1861
          CALL FLUSH(6)
1862
        END IF
1863

1864
        IF( Robust ) THEN
1865
          ! Always store the best solution so far (with some small margin)
1866
          IF( errorInd < RobustStep * BestNorm ) THEN
1867
            BestIter = iter
1868
            BestNorm = errorInd
1869
            IF (Smoothing) THEN
1870
              Bestx = x_s
1871
            ELSE
1872
              Bestx = x
1873
            END IF
1874
            BadIterCount = 0
1875
          ELSE
1876
            BadIterCount = BadIterCount + 1
1877
          END IF
1878

1879
          ! If we have diverged too much and have found already a good candidate, then take it
1880
          IF( BestNorm <  RobustTol .AND. &
1881
              ( errorInd > RobustMaxTol .OR. BadIterCount > MaxBadIter ) ) THEN
1882
            EXIT
1883
          END IF
1884
          
1885
        END IF
1886
                        
1887
        Converged = (errorind < Tol)
1888
        Diverged = (errorind > MaxTol) .OR. (errorind /= errorind)
1889
        IF (iter == MaxRounds) EXIT
1890
      END DO ! end of while loop
1891

1892
      IF( Smoothing ) x = x_s
1893
      
1894
      IF( Robust ) THEN
1895
        IF( BestNorm < RobustTol ) THEN
1896
          Converged = .TRUE.
1897
        END IF
1898
        IF( BestNorm < errorInd ) THEN
1899
          x = Bestx
1900
        END IF        
1901
        IF(OutputInterval /= HUGE(OutputInterval)) THEN
1902
          WRITE(*,'(A,I8,E11.4,I8,E11.4)') 'Idrs robust: ',&
1903
              iter, BestNorm, BestIter, errorind
1904
          CALL FLUSH(6)
1905
        END IF
1906
      ELSE
1907
        IF(OutputInterval /= HUGE(OutputInterval)) THEN
1908
          WRITE(*,'(A,I8,E11.4)') 'Idrs: ',&
1909
              iter, errorind
1910
          CALL FLUSH(6)
1911
        END IF
1912
      END IF
1913
      
1914
    !----------------------------------------------------------
1915
    END SUBROUTINE RealIDRS
1916
    !----------------------------------------------------------
1917

1918
!--------------------------------------------------------------
1919
  END SUBROUTINE itermethod_idrs
1920
!--------------------------------------------------------------
1921

1922

1923
!------------------------------------------------------------------------------
1924
!> This routine provides the complex version to the GCR linear solver.
1925
!------------------------------------------------------------------------------
1926
 SUBROUTINE itermethod_z_gcr( xvec, rhsvec, &
1927
      ipar, dpar, work, matvecsubr, pcondlsubr, &
1928
      pcondrsubr, dotprodfun, normfun, stopcfun )
1929
!------------------------------------------------------------------------------
1930
    USE huti_interfaces
1931
    IMPLICIT NONE
1932
    PROCEDURE( mv_iface_z ), POINTER :: matvecsubr
1933
    PROCEDURE( pc_iface_z ), POINTER :: pcondlsubr
1934
    PROCEDURE( pc_iface_z ), POINTER :: pcondrsubr
1935
    PROCEDURE( dotp_iface_z ), POINTER :: dotprodfun
1936
    PROCEDURE( norm_iface_z ), POINTER :: normfun
1937
    PROCEDURE( stopc_iface_z ), POINTER :: stopcfun
1938

1939
    INTEGER, DIMENSION(HUTI_IPAR_DFLTSIZE) :: ipar
1940
    COMPLEX(KIND=dp), TARGET, DIMENSION(HUTI_NDIM) :: xvec, rhsvec
1941
    DOUBLE PRECISION, DIMENSION(HUTI_DPAR_DFLTSIZE) :: dpar
1942
    COMPLEX(KIND=dp), DIMENSION(HUTI_WRKDIM,HUTI_NDIM) :: work
1943

1944
    COMPLEX(KIND=dp) :: y(HUTI_NDIM),f(HUTI_NDIM)
1945
    INTEGER :: ndim, RestartN, i
1946
    INTEGER :: Rounds, OutputInterval
1947
    REAL(KIND=dp) :: MinTol, MaxTol, Residual
1948
    LOGICAL :: Converged, Diverged, UseStopCFun
1949

1950
    CALL Info('Itermethod_z_gcr','Starting GCR iteration',Level=25)
1951
    
1952
    ndim = HUTI_NDIM
1953
    Rounds = HUTI_MAXIT
1954
    MinTol = HUTI_TOLERANCE
1955
    MaxTol = HUTI_MAXTOLERANCE
1956
    OutputInterval = HUTI_DBUGLVL
1957
    RestartN = HUTI_GCR_RESTART 
1958
    UseStopCFun = HUTI_STOPC == HUTI_USUPPLIED_STOPC
1959
    
1960
    Converged = .FALSE.
1961
    Diverged = .FALSE.
1962
    
1963
    !----------------------------------------------------------------------------
1964
    ! Transform the solution vector and the right-hand side vector to 
1965
    ! complex-valued vectors y and f
1966
    !---------------------------------------------------------------------------
1967
    DO i=1,ndim
1968
      y(i)=xvec(i)
1969
      f(i)=rhsvec(i)
1970
    END DO
1971
       
1972
    CALL GCR_Z(ndim, GlobalMatrix, y, f, Rounds, MinTol, MaxTol, Residual, &
1973
        Converged, Diverged, OutputInterval, RestartN )
1974

1975
    IF(Converged) HUTI_INFO = HUTI_CONVERGENCE
1976
    IF(Diverged) HUTI_INFO = HUTI_DIVERGENCE
1977
    IF ( (.NOT. Converged) .AND. (.NOT. Diverged) ) HUTI_INFO = HUTI_MAXITER
1978
   
1979
    DO i=1,ndim
1980
      xvec(i) = y(i)
1981
    END DO
1982

1983
  CONTAINS 
1984
    
1985
    
1986
!------------------------------------------------------------------------------  
1987
    SUBROUTINE GCR_Z( n, A, x, b, Rounds, MinTolerance, MaxTolerance, Residual, &
1988
        Converged, Diverged, OutputInterval, m) 
1989
!------------------------------------------------------------------------------
1990
      TYPE(Matrix_t), POINTER :: A
1991
      INTEGER :: Rounds
1992
      COMPLEX(KIND=dp) :: x(n),b(n)
1993
      LOGICAL :: Converged, Diverged
1994
      REAL(KIND=dp) :: MinTolerance, MaxTolerance, Residual
1995
      INTEGER :: n, OutputInterval, m
1996
!------------------------------------------------------------------------------      
1997
      REAL(KIND=dp) :: bnorm,rnorm
1998
      COMPLEX(KIND=dp), ALLOCATABLE :: R(:)
1999
      COMPLEX(KIND=dp), ALLOCATABLE :: S(:,:), V(:,:), T1(:), T2(:)
2000
      INTEGER :: i,j,k,allocstat
2001
      COMPLEX(KIND=dp) :: beta, czero
2002
      REAL(KIND=dp) :: alpha, trueresnorm, normerr
2003
      COMPLEX(KIND=dp), ALLOCATABLE :: trueres(:)
2004
!------------------------------------------------------------------------------
2005
            
2006
      ALLOCATE( R(n), T1(n), T2(n), trueres(n), STAT=allocstat )
2007
      IF( allocstat /= 0) &
2008
          CALL Fatal('GCR_Z','Failed to allocate memory of size: '//I2S(n))
2009
      IF ( m > 1 ) THEN        
2010
         ALLOCATE( S(n,m-1), V(n,m-1), STAT=allocstat )
2011
         IF ( allocstat /= 0 ) THEN
2012
           CALL Fatal('GCR_Z','Failed to allocate memory of size: '&
2013
               //I2S(n)//' x '//I2S(m-1))
2014
         END IF
2015
         czero = CMPLX( 0.0_dp, 0.0_dp )
2016
         V(1:n,1:m-1) = czero
2017
         S(1:n,1:m-1) = czero
2018
      END IF	
2019
      
2020
      CALL matvecsubr( x, r, ipar )
2021
      r(1:n) = b(1:n) - r(1:n)
2022
      
2023
      bnorm = normfun(n, b, 1)
2024
      rnorm = normfun(n, r, 1)
2025

2026
      IF (UseStopCFun) THEN
2027
        Residual = stopcfun(x,b,r,ipar,dpar)
2028
      ELSE
2029
        Residual = rnorm / bnorm
2030
      END IF
2031
      Converged = (Residual < MinTolerance) 
2032
      Diverged = (Residual > MaxTolerance) .OR. (Residual /= Residual)    
2033
      IF( Converged .OR. Diverged) RETURN
2034
      
2035
      DO k=1,Rounds
2036
	 !----------------------------------------------
2037
	 ! Check for restarting
2038
         !--------------------------------------------- 
2039
         IF ( MOD(k,m)==0 ) THEN
2040
            j = m
2041
         ELSE
2042
            j = MOD(k,m)
2043
            !--------------------------------------------
2044
            ! Compute the true residual when restarting:
2045
            !--------------------------------------------
2046
            IF ( (j==1) .AND. (k>1) ) THEN
2047
               CALL matvecsubr( x, r, ipar ) 
2048
               r(1:n) = b(1:n) - r(1:n)
2049
            END IF
2050
         END IF
2051
         !----------------------------------------------------------
2052
         ! Perform the preconditioning...
2053
         !---------------------------------------------------------------
2054
         CALL pcondrsubr( T1, r, ipar )         
2055
         CALL matvecsubr( T1, T2, ipar )
2056
         !--------------------------------------------------------------
2057
         ! Perform the orthogonalization of the search directions....
2058
         !--------------------------------------------------------------
2059
         DO i=1,j-1
2060
            beta = dotprodfun(n, V(1:n,i), 1, T2(1:n), 1 )
2061

2062
            T1(1:n) = T1(1:n) - beta * S(1:n,i)
2063
            T2(1:n) = T2(1:n) - beta * V(1:n,i)        
2064
         END DO
2065

2066
         alpha = normfun(n, T2(1:n), 1 )
2067
         T1(1:n) = 1.0d0/alpha * T1(1:n)
2068
         T2(1:n) = 1.0d0/alpha * T2(1:n)
2069

2070
         !-------------------------------------------------------------
2071
         ! The update of the solution and save the search data...
2072
         !------------------------------------------------------------- 
2073
         beta = dotprodfun(n, T2(1:n), 1, r(1:n), 1 )
2074
         x(1:n) = x(1:n) + beta * T1(1:n)      
2075
         r(1:n) = r(1:n) - beta * T2(1:n)
2076
	 IF ( j /= m ) THEN
2077
            S(1:n,j) = T1(1:n)
2078
            V(1:n,j) = T2(1:n)
2079
	 END IF       
2080

2081
         !--------------------------------------------------------------
2082
         ! Check whether the convergence criterion is met 
2083
         !--------------------------------------------------------------
2084
         rnorm = normfun(n, r, 1)
2085

2086
         IF (UseStopCFun) THEN
2087
           Residual = stopcfun(x,b,r,ipar,dpar)
2088
           IF( MOD(k,OutputInterval) == 0) THEN
2089
             WRITE (*, '(A, I6, 2E12.4)') '   gcr:',k, rnorm / bnorm, residual
2090
             CALL FLUSH(6)
2091
           END IF
2092
         ELSE
2093
           Residual = rnorm / bnorm
2094
           IF( MOD(k,OutputInterval) == 0) THEN
2095
             WRITE (*, '(A, I8, 3ES12.4,A)') '   gcrz:',k, residual, beta,'i'
2096
             CALL FLUSH(6)
2097
           END IF
2098
         END IF
2099
        
2100
         Converged = (Residual < MinTolerance)
2101
         !-----------------------------------------------------------------
2102
         ! Make an additional check that the true residual agrees with 
2103
         ! the iterated residual:
2104
         !-----------------------------------------------------------------
2105
         IF (Converged ) THEN
2106
            CALL matvecsubr( x, trueres, ipar )
2107
            trueres(1:n) = b(1:n) - trueres(1:n)
2108
            TrueResNorm = normfun(n, trueres, 1)
2109
            NormErr = ABS(TrueResNorm - rnorm)/TrueResNorm
2110
            IF ( NormErr > 1.0d-1 ) THEN
2111
               CALL Info('WARNING', 'Iterated GCR solution may not be accurate', Level=2)
2112
               WRITE( Message, * ) 'Iterated GCR residual norm = ', rnorm
2113
               CALL Info('WARNING', Message, Level=2)
2114
               WRITE( Message, * ) 'True residual norm = ', TrueResNorm
2115
               CALL Info('WARNING', Message, Level=2)   
2116
               CALL FLUSH(6)
2117
             END IF
2118
         END IF 
2119
         Diverged = (Residual > MaxTolerance) .OR. (Residual /= Residual)    
2120
         IF( Converged .OR. Diverged) EXIT
2121
        
2122
      END DO
2123
      
2124
      DEALLOCATE( R, T1, T2 )
2125
      IF ( m > 1 ) DEALLOCATE( S, V)
2126
      
2127
    END SUBROUTINE GCR_Z
2128
    
2129
!------------------------------------------------------------------------------
2130
  END SUBROUTINE itermethod_z_gcr
2131
!------------------------------------------------------------------------------
2132

2133

2134

2135
!------------------------------------------------------------------------------
2136
!> This routine provides a complex version of the BiCGstab(l) solver for linear systems.
2137
!------------------------------------------------------------------------------
2138
  SUBROUTINE itermethod_z_bicgstabl( xvec, rhsvec, &
2139
       ipar, dpar, work, matvecsubr, pcondlsubr, &
2140
       pcondrsubr, dotprodfun, normfun, stopcfun )
2141
 !------------------------------------------------------------------------------
2142
    USE huti_interfaces
2143
    IMPLICIT NONE
2144
    PROCEDURE( mv_iface_z ), POINTER :: matvecsubr
2145
    PROCEDURE( pc_iface_z ), POINTER :: pcondlsubr
2146
    PROCEDURE( pc_iface_z ), POINTER :: pcondrsubr
2147
    PROCEDURE( dotp_iface_z ), POINTER :: dotprodfun
2148
    PROCEDURE( norm_iface_z ), POINTER :: normfun
2149
    PROCEDURE( stopc_iface_z ), POINTER :: stopcfun
2150

2151
    INTEGER, DIMENSION(HUTI_IPAR_DFLTSIZE) :: ipar
2152
    COMPLEX(KIND=dp), TARGET, DIMENSION(HUTI_NDIM) :: xvec, rhsvec
2153
    DOUBLE PRECISION, DIMENSION(HUTI_DPAR_DFLTSIZE) :: dpar
2154
    COMPLEX(KIND=dp), DIMENSION(HUTI_WRKDIM,HUTI_NDIM) :: work
2155

2156
    COMPLEX(KIND=dp) :: y(HUTI_NDIM),f(HUTI_NDIM)
2157
    INTEGER :: ndim, i, PolynomialDegree
2158
    INTEGER :: Rounds, OutputInterval
2159
    REAL(KIND=dp) :: MinTol, MaxTol
2160
    LOGICAL :: Converged, Diverged
2161
    !--------------------------------------------------------------------------------    
2162

2163
    ndim = HUTI_NDIM
2164
    Rounds = HUTI_MAXIT
2165
    MinTol = HUTI_TOLERANCE
2166
    MaxTol = HUTI_MAXTOLERANCE
2167
    OutputInterval = HUTI_DBUGLVL
2168
    PolynomialDegree = HUTI_BICGSTABL_L 
2169
    !----------------------------------------------------------------------------
2170
    ! Transform the solution vector and the right-hand side vector to 
2171
    ! complex-valued vectors y and f
2172
    !---------------------------------------------------------------------------
2173
    DO i=1,ndim
2174
      y(i)=xvec(i)
2175
      f(i)=rhsvec(i)
2176
    END DO
2177

2178
    CALL ComplexBiCGStabl(ndim, GlobalMatrix, y, f, Rounds, MinTol, MaxTol, &
2179
         Converged, Diverged, OutputInterval, PolynomialDegree)
2180

2181
    IF(Converged) HUTI_INFO = HUTI_CONVERGENCE
2182
    IF(Diverged) HUTI_INFO = HUTI_DIVERGENCE
2183
    IF ( (.NOT. Converged) .AND. (.NOT. Diverged) ) HUTI_INFO = HUTI_MAXITER
2184

2185
    DO i=1,ndim
2186
      xvec(i) = y(i)
2187
    END DO
2188

2189
  CONTAINS 
2190

2191
    !-----------------------------------------------------------------------------------
2192
    SUBROUTINE ComplexBiCGStabl( n, A, x, b, MaxRounds, Tol, MaxTol, Converged, &
2193
         Diverged, OutputInterval, l)
2194
    !-----------------------------------------------------------------------------------
2195
    !   This subroutine solves complex linear systems Ax = b by using the BiCGStab(l) algorithm 
2196
    !   with l >= 2 and the right-oriented preconditioning. 
2197
    !
2198
    !   The subroutine has been written using as a starting point the work of D.R. Fokkema 
2199
    !   (subroutine zbistbl v1.1 1998). Dr. Fokkema has given the right to distribute
2200
    !   the derived work under GPL and hence the original copyright notice of the subroutine
2201
    !   has been removed accordingly.  
2202
    !
2203
    !----------------------------------------------------------------------------------- 
2204
      INTEGER :: l   ! polynomial degree
2205
      INTEGER :: n, MaxRounds, OutputInterval   
2206
      LOGICAL :: Converged, Diverged
2207
      TYPE(Matrix_t), POINTER :: A
2208
      COMPLEX(KIND=dp) :: x(n), b(n)
2209
      REAL(KIND=dp) :: Tol, MaxTol
2210
      !------------------------------------------------------------------------------
2211
      COMPLEX(KIND=dp) :: zzero, zone, t(n), kappa0, kappal 
2212
      REAL(KIND=dp) :: rnrm0, rnrm, mxnrmx, mxnrmr, errorind, &
2213
           delta = 1.0d-2, bnrm
2214
      INTEGER :: i, j, rr, r, u, xp, bp, z, zz, y0, yl, y, k, iwork(l-1), stat, Round
2215
      COMPLEX(KIND=dp) :: alpha, beta, omega, rho0, rho1, sigma, zdotc, varrho, hatgamma
2216
      COMPLEX(KIND=dp), ALLOCATABLE :: work(:,:), rwork(:,:)
2217
      LOGICAL rcmp, xpdt, EarlyExit
2218
      !------------------------------------------------------------------------------
2219

2220
      IF ( l < 2) CALL Fatal( 'RealBiCGStabl', 'Polynomial degree < 2' )
2221

2222
      IF ( ALL(x == CMPLX(0.0d0,0.0d0,kind=dp)) ) x = b
2223

2224
      zzero = CMPLX( 0.0d0,0.0d0, KIND=dp)
2225
      zone =  CMPLX( 1.0d0,0.0d0, KIND=dp)
2226

2227
      ALLOCATE( work(n,3+2*(l+1)), rwork(l+1,3+2*(l+1)) )
2228
      work = CMPLX( 0.0d0, 0.0d0, KIND=dp )
2229
      rwork = CMPLX( 0.0d0, 0.0d0, KIND=dp )
2230

2231
      rr = 1
2232
      r = rr+1
2233
      u = r+(l+1)
2234
      xp = u+(l+1)
2235
      bp = xp+1
2236

2237
      z = 1
2238
      zz = z+(l+1)
2239
      y0 = zz+(l+1)
2240
      yl = y0+1
2241
      y = yl+1
2242

2243
      CALL matvecsubr( x, work(1:n,r), ipar )
2244
      work(1:n,r) = b(1:n) - work(1:n,r)
2245
      bnrm = normfun(n, b(1:n), 1)
2246
      rnrm0 = normfun(n, work(1:n,r), 1)
2247

2248
      !-------------------------------------------------------------------
2249
      ! Check whether the initial guess satisfies the stopping criterion
2250
      !--------------------------------------------------------------------
2251
      errorind = rnrm0 / bnrm
2252
      Converged = (errorind < Tol)
2253
      Diverged = (errorind > MaxTol) .OR. (errorind /= errorind)
2254

2255
      IF( Converged .OR. Diverged) RETURN
2256
      EarlyExit = .FALSE.
2257

2258
      work(1:n,rr) = work(1:n,r) 
2259
      work(1:n,bp) = work(1:n,r)
2260
      work(1:n,xp) = x(1:n)
2261

2262
      rnrm = rnrm0
2263
      mxnrmx = rnrm0
2264
      mxnrmr = rnrm0
2265
      x(1:n) = zzero    
2266
      alpha = zzero
2267
      omega = zone
2268
      sigma = zone
2269
      rho0 = zone
2270

2271
      DO Round=1,MaxRounds 
2272
         !-------------------------
2273
         ! --- The BiCG part ---
2274
         !-------------------------
2275
         rho0 = -omega*rho0
2276

2277
         DO k=1,l
2278
            rho1 = dotprodfun(n, work(1:n,rr), 1, work(1:n,r+k-1), 1)
2279
            IF (rho0 == zzero) THEN
2280
               CALL Fatal( 'ComplexBiCGStab(l)', 'Breakdown error 1.' )
2281
            ENDIF
2282
            beta = alpha*(rho1/rho0)
2283
            rho0 = rho1
2284
            DO j=0,k-1
2285
               work(1:n,u+j) = work(1:n,r+j) - beta*work(1:n,u+j)
2286
            ENDDO
2287
            CALL pcondrsubr( t, work(1:n,u+k-1), ipar )
2288
            CALL matvecsubr( t, work(1:n,u+k),   ipar )
2289

2290
            sigma = dotprodfun(n, work(1:n,rr), 1, work(1:n,u+k), 1)
2291
            IF (sigma == zzero) THEN
2292
               CALL Fatal( 'ComplexBiCGStab(l)', 'Breakdown error 2.' )
2293
            ENDIF
2294
            alpha = rho1/sigma
2295
            x(1:n) = x(1:n) + alpha * work(1:n,u)
2296
            DO j=0,k-1
2297
               work(1:n,r+j) = work(1:n,r+j) - alpha * work(1:n,u+j+1)
2298
            ENDDO
2299
            CALL pcondrsubr( t, work(1:n,r+k-1), ipar )
2300
            CALL matvecsubr( t, work(1:n,r+k),   ipar )
2301
            rnrm = normfun(n, work(1:n,r), 1)
2302
            mxnrmx = MAX (mxnrmx, rnrm)
2303
            mxnrmr = MAX (mxnrmr, rnrm)
2304

2305
            !----------------------------------------------------------------------
2306
            ! In some simple cases, a few BiCG updates may already be enough to
2307
            ! obtain the solution. The following is for handling this special case. 
2308
            !----------------------------------------------------------------------
2309
            errorind = rnrm / bnrm
2310
            Converged = (errorind < Tol) 
2311
            IF (Converged) THEN
2312
               EarlyExit = .TRUE.
2313
               EXIT
2314
            END IF
2315

2316
         ENDDO
2317
         
2318
         IF (EarlyExit) EXIT        
2319

2320
         !--------------------------------------
2321
         ! --- The convex polynomial part ---
2322
         !--------------------------------------
2323

2324
         DO i=1,l+1
2325
            DO j=1,i
2326
               rwork(i,j) = dotprodfun(n, work(1:n,r+i-1), 1, work(1:n,r+j-1),1 ) 
2327
            END DO
2328
         END DO
2329
         DO j=2,l+1
2330
            rwork(1:j-1,j) = CONJG( rwork(j,1:j-1) )
2331
         END DO
2332

2333
         rwork(1:l+1,zz:zz+l) = rwork(1:l+1,z:z+l)
2334
         CALL zgetrf (l-1, l-1, rwork(2:l,zz+1:zz+l-1), l-1, &
2335
              iwork, stat)
2336

2337
         ! --- tilde r0 and tilde rl (small vectors)
2338

2339
         rwork(1,y0) = -zone
2340
         rwork(2:l,y0) = rwork(2:l,z) 
2341
         CALL zgetrs('n', l-1, 1, rwork(2:l,zz+1:zz+l-1), l-1, iwork, &
2342
              rwork(2:l,y0), l-1, stat)
2343
         rwork(l+1,y0) = zzero
2344

2345
         rwork(1,yl) = zzero
2346
         rwork(2:l,yl) = rwork(2:l,z+l) 
2347
         CALL zgetrs ('n', l-1, 1, rwork(2:l,zz+1:zz+l-1), l-1, iwork, &
2348
              rwork(2:l,yl), l-1, stat)
2349
         rwork(l+1,yl) = -zone
2350

2351
         ! --- Convex combination
2352
         call zmv( rwork(1:l+1,z:z+l), rwork(1:l+1,y0), rwork(1:l+1,y), l+1 )
2353
         kappa0 = SQRT( ABS(zdotc(l+1, rwork(1:l+1,y0), 1, rwork(1:l+1,y), 1)) )  ! replace zdotc
2354

2355
         call zmv( rwork(1:l+1,z:z+l), rwork(1:l+1,yl), rwork(1:l+1,y), l+1 )
2356
         kappal = SQRT( ABS(zdotc(l+1, rwork(1:l+1,yl), 1, rwork(1:l+1,y), 1)) )  ! replace zdotc
2357

2358
         call zmv( rwork(1:l+1,z:z+l), rwork(1:l+1,y0), rwork(1:l+1,y), l+1 )
2359
         varrho = zdotc(l+1, rwork(1:l+1,yl), 1, rwork(1:l+1,y), 1) / &           ! replace zdotc
2360
              (kappa0*kappal)
2361

2362
         hatgamma = varrho/ABS(varrho) * MAX(ABS(varrho),7d-1) * &
2363
              kappa0/kappal
2364
         rwork(1:l+1,y0) = rwork(1:l+1,y0) - hatgamma * rwork(1:l+1,yl)
2365

2366
         !  --- Update
2367

2368
         omega = rwork(l+1,y0)
2369
         DO j=1,l
2370
            work(1:n,u) = work(1:n,u) - rwork(j+1,y0) * work(1:n,u+j)
2371
            x(1:n) = x(1:n) + rwork(j+1,y0) * work(1:n,r+j-1)
2372
            work(1:n,r) = work(1:n,r) - rwork(j+1,y0) * work(1:n,r+j)
2373
         ENDDO
2374

2375
         call zmv( rwork(1:l+1,z:z+l), rwork(1:l+1,y0), rwork(1:l+1,y), l+1 )
2376
         rnrm = SQRT( ABS(zdotc(l+1, rwork(1:l+1,y0), 1, rwork(1:l+1,y), 1)) )
2377

2378
         !---------------------------------------
2379
         !  --- The reliable update part ---
2380
         !---------------------------------------
2381

2382
         mxnrmx = MAX (mxnrmx, rnrm)
2383
         mxnrmr = MAX (mxnrmr, rnrm)
2384
         xpdt = (rnrm < delta*rnrm0 .AND. rnrm0 < mxnrmx)
2385
         rcmp = ((rnrm < delta*mxnrmr .AND. rnrm0 < mxnrmr) .OR. xpdt)
2386
         IF (rcmp) THEN
2387
            ! PRINT *, 'Performing residual update...'
2388
            CALL pcondrsubr( t, x, ipar )
2389
            CALL matvecsubr( t, work(1:n,r), ipar )
2390
            work(1:n,r) = work(1:n,bp) - work(1:n,r)
2391
            mxnrmr = rnrm
2392
            IF (xpdt) THEN
2393
               ! PRINT *, 'Performing solution update...'
2394
               work(1:n,xp) = work(1:n,xp) + t(1:n)
2395
               x(1:n) = zzero
2396
               work(1:n,bp) = work(1:n,r)
2397
               mxnrmx = rnrm
2398
            ENDIF
2399
         ENDIF
2400

2401
         IF (rcmp) THEN
2402
            IF (xpdt) THEN       
2403
               t(1:n) = work(1:n,xp)
2404
            ELSE
2405
               t(1:n) = t(1:n) + work(1:n,xp)  
2406
            END IF
2407
         ELSE
2408
            CALL pcondrsubr( t, x, ipar )
2409
            t(1:n) = t(1:n)+work(1:n,xp)
2410
         END IF
2411

2412
         errorind = rnrm/bnrm
2413
         IF( MOD(Round,OutputInterval) == 0) THEN
2414
           WRITE (*, '(I8, E11.4)') Round, errorind
2415
           CALL FLUSH(6)
2416
         END IF
2417

2418
         Converged = (errorind < Tol) 
2419
         Diverged = (errorind > MaxTol) .OR. (errorind /= errorind)    
2420
         IF( Converged .OR. Diverged) EXIT    
2421
      END DO
2422

2423
      IF( EarlyExit .AND. (OutputInterval/=HUGE(OutputInterval)) ) THEN
2424
        WRITE (*, '(I8, E11.4)') Round, errorind
2425
        CALL FLUSH(6)
2426
      END IF
2427

2428
      !------------------------------------------------------------
2429
      ! We have solved z = P*x, so finally solve the true unknown x
2430
      !------------------------------------------------------------
2431
      t(1:n) = x(1:n)
2432
      CALL pcondrsubr( x, t, ipar )
2433
      x(1:n) = x(1:n) + work(1:n,xp)      
2434

2435
    !----------------------------------------------------------
2436
    END SUBROUTINE ComplexBiCGStabl
2437
    !----------------------------------------------------------
2438

2439
!--------------------------------------------------------------
2440

2441
    SUBROUTINE zmv(a,x,y,n)
2442
      integer, INTENT(in) :: n
2443
      complex(kind=dp), INTENT(in)  ::a(n,n), x(n)
2444
      complex(kind=dp), INTENT(out) ::y(n)
2445
 
2446
      complex(kind=dp), parameter :: zone  = cmplx(1._dp, 0._dp,kind=dp)
2447
      complex(kind=dp), parameter :: zzero = cmplx(0._dp, 0._dp,kind=dp)
2448
    
2449
      IF(n>8) THEN
2450
        CALL zhemv ('u', n, zone, a, n, x, 1, zzero, y, 1)
2451
        RETURN
2452
      END IF
2453

2454
      y(1) = a(1,1)*x(1) + a(1,2)*x(2) + a(1,3)*x(3)
2455
      y(2) = a(2,1)*x(1) + a(2,2)*x(2) + a(2,3)*x(3)
2456
      y(3) = a(3,1)*x(1) + a(3,2)*x(2) + a(3,3)*x(3)
2457
      IF(n==3) RETURN
2458

2459
      y(1) = y(1) + a(1,4)*x(4)
2460
      y(2) = y(2) + a(2,4)*x(4)
2461
      y(3) = y(3) + a(3,4)*x(4)
2462
      y(4) = a(4,1)*x(1)+a(4,2)*x(2)+a(4,3)*x(3)+a(4,4)*x(4)
2463
      IF(n==4) RETURN
2464

2465
      y(1) = y(1) + a(1,5)*x(5)
2466
      y(2) = y(2) + a(2,5)*x(5)
2467
      y(3) = y(3) + a(3,5)*x(5)
2468
      y(4) = y(4) + a(4,5)*x(5)
2469
      y(5) = a(5,1)*x(1)+a(5,2)*x(2)+a(5,3)*x(3)+a(5,4)*x(4)+ &
2470
             a(5,5)*x(5)
2471
      IF(n==5) RETURN
2472

2473
      y(1) = y(1) + a(1,6)*x(6)
2474
      y(2) = y(2) + a(2,6)*x(6)
2475
      y(3) = y(3) + a(3,6)*x(6)
2476
      y(4) = y(4) + a(4,6)*x(6)
2477
      y(5) = y(5) + a(5,6)*x(6)
2478
      y(6) = a(6,1)*x(1)+a(6,2)*x(2)+a(6,3)*x(3)+a(6,4)*x(4)+ &
2479
             a(6,5)*x(5)+a(6,6)*x(6)
2480
      IF(n==6) RETURN
2481

2482
      y(1) = y(1) + a(1,7)*x(7)
2483
      y(2) = y(2) + a(2,7)*x(7)
2484
      y(3) = y(3) + a(3,7)*x(7)
2485
      y(4) = y(4) + a(4,7)*x(7)
2486
      y(5) = y(5) + a(5,7)*x(7)
2487
      y(6) = y(6) + a(6,7)*x(7)
2488
      y(7) = a(7,1)*x(1)+a(7,2)*x(2)+a(7,3)*x(3)+a(7,4)*x(4)+ &
2489
             a(7,5)*x(5)+a(7,6)*x(6)+a(7,7)*x(7)
2490
      IF(n==7) RETURN
2491

2492
      y(1) = y(1) + a(1,8)*x(8)
2493
      y(2) = y(2) + a(2,8)*x(8)
2494
      y(3) = y(3) + a(3,8)*x(8)
2495
      y(4) = y(4) + a(4,8)*x(8)
2496
      y(5) = y(5) + a(5,8)*x(8)
2497
      y(6) = y(6) + a(6,8)*x(8)
2498
      y(7) = y(7) + a(7,8)*x(8)
2499
      y(8) = a(8,1)*x(1)+a(8,2)*x(2)+a(8,3)*x(3)+a(8,4)*x(4)+ &
2500
             a(8,5)*x(5)+a(8,6)*x(6)+a(8,7)*x(7)+a(8,8)*x(8)
2501

2502
    END SUBROUTINE zmv
2503

2504
!--------------------------------------------------------------
2505
  END SUBROUTINE itermethod_z_bicgstabl
2506
!--------------------------------------------------------------
2507

2508

2509
!------------------------------------------------------------------------------
2510
!> This routine provides a complex version of the IDR(s) solver for linear systems.
2511
!------------------------------------------------------------------------------
2512
  SUBROUTINE itermethod_z_idrs( xvec, rhsvec, &
2513
      ipar, dpar, work, matvecsubr, pcondlsubr, &
2514
      pcondrsubr, dotprodfun, normfun, stopcfun )
2515
!------------------------------------------------------------------------------
2516
    USE huti_interfaces
2517
    IMPLICIT NONE
2518
    PROCEDURE( mv_iface_z ), POINTER :: matvecsubr
2519
    PROCEDURE( pc_iface_z ), POINTER :: pcondlsubr
2520
    PROCEDURE( pc_iface_z ), POINTER :: pcondrsubr
2521
    PROCEDURE( dotp_iface_z ), POINTER :: dotprodfun
2522
    PROCEDURE( norm_iface_z ), POINTER :: normfun
2523
    PROCEDURE( stopc_iface_z ), POINTER :: stopcfun
2524

2525
    INTEGER, DIMENSION(HUTI_IPAR_DFLTSIZE) :: ipar
2526
    COMPLEX(KIND=dp), TARGET, DIMENSION(HUTI_NDIM) :: xvec, rhsvec
2527
    DOUBLE PRECISION, DIMENSION(HUTI_DPAR_DFLTSIZE) :: dpar
2528
    COMPLEX(KIND=dp), DIMENSION(HUTI_WRKDIM,HUTI_NDIM) :: work
2529

2530
    COMPLEX(KIND=dp) :: y(HUTI_NDIM),f(HUTI_NDIM)
2531
    INTEGER :: ndim, i, s
2532
    INTEGER :: Rounds, OutputInterval
2533
    REAL(KIND=dp) :: MinTol, MaxTol
2534
    LOGICAL :: Converged, Diverged
2535
    !--------------------------------------------------------------------------------    
2536

2537
    ndim = HUTI_NDIM
2538
    Rounds = HUTI_MAXIT
2539
    MinTol = HUTI_TOLERANCE
2540
    MaxTol = HUTI_MAXTOLERANCE
2541
    OutputInterval = HUTI_DBUGLVL
2542
    s = HUTI_IDRS_S 
2543
    !----------------------------------------------------------------------------
2544
    ! Transform the solution vector and the right-hand side vector to 
2545
    ! complex-valued vectors y and f
2546
    !---------------------------------------------------------------------------
2547
    DO i=1,ndim
2548
      y(i)=xvec(i)
2549
      f(i)=rhsvec(i)
2550
    END DO
2551
    CALL ComplexIDRS(ndim, GlobalMatrix, y, f, Rounds, MinTol, MaxTol, &
2552
         Converged, Diverged, OutputInterval, s )
2553

2554
    IF(Converged) HUTI_INFO = HUTI_CONVERGENCE
2555
    IF(Diverged) HUTI_INFO = HUTI_DIVERGENCE
2556
    IF ( (.NOT. Converged) .AND. (.NOT. Diverged) ) HUTI_INFO = HUTI_MAXITER
2557

2558
    DO i=1,ndim
2559
      xvec(i) = y(i)
2560
    END DO
2561

2562
  CONTAINS 
2563

2564
!----------------------------------------------------------------------------------- 
2565
!   This subroutine solves complex linear systems Ax = b by using the IDR(s) algorithm 
2566
!   with s >= 1 and the right-oriented preconditioning. 
2567
!
2568
!   The subroutine ComplexIDRS has been written by M.B. van Gijzen
2569
!-----------------------------------------------------------------------------------
2570
    SUBROUTINE ComplexIDRS( n, A, x, b, MaxRounds, Tol, MaxTol, Converged, &
2571
        Diverged, OutputInterval, s )
2572
!-----------------------------------------------------------------------------------
2573
      INTEGER :: s  
2574
      INTEGER :: n, MaxRounds, OutputInterval   
2575
      LOGICAL :: Converged, Diverged, UseStopCFun
2576
      TYPE(Matrix_t), POINTER :: A
2577
      COMPLEX(KIND=dp) :: x(n), b(n)
2578
      REAL(KIND=dp) :: Tol, MaxTol
2579
!------------------------------------------------------------------------------
2580

2581
      ! Local arrays:
2582
      REAL(kind=dp) :: Pr(n,s), Pi(n,s) 
2583
      COMPLEX(kind=dp) :: P(n,s)
2584
      COMPLEX(kind=dp) :: G(n,s)
2585
      COMPLEX(kind=dp) :: U(n,s)
2586
      COMPLEX(kind=dp) :: r(n) 
2587
      COMPLEX(kind=dp) :: v(n)   
2588
      COMPLEX(kind=dp) :: t(n)  
2589
      COMPLEX(kind=dp) :: M(s,s), f(s), mu(s)
2590
      COMPLEX(kind=dp) :: alpha(s), beta(s), gamma(s)
2591

2592
      COMPLEX(kind=dp) :: om, tr    
2593
      REAL(kind=dp) :: nr, nt, rho, kappa
2594

2595
      INTEGER :: iter                         ! number of iterations
2596
      INTEGER :: ii                           ! inner iterations index
2597
      INTEGER :: jj                           ! G-space index
2598
      REAL(kind=dp) :: normb, normr, errorind ! for tolerance check
2599
      INTEGER :: i,j,k,l                      ! loop counters
2600

2601
      UseStopCFun = HUTI_STOPC == HUTI_USUPPLIED_STOPC
2602
      
2603
      U = 0.0d0
2604

2605
      ! Compute initial residual, set absolute tolerance
2606
      normb = normfun(n,b,1)
2607
      CALL matvecsubr( x, t, ipar )
2608
      r = b - t
2609
      IF (UseStopCFun) THEN
2610
        errorind = stopcfun(x,b,r,ipar,dpar)
2611
      ELSE
2612
        normr = normfun(n,r,1)
2613
        errorind = normr / normb
2614
      END IF
2615
      
2616
      !-------------------------------------------------------------------
2617
      ! Check whether the initial guess satisfies the stopping criterion
2618
      !--------------------------------------------------------------------
2619
      Converged = (errorind < Tol)
2620
      Diverged = (errorind > MaxTol) .OR. (errorind /= errorind)
2621

2622
      IF ( Converged .OR. Diverged ) RETURN
2623

2624
      ! Define P and kappa 
2625
      CALL RANDOM_NUMBER(Pr)
2626
      CALL RANDOM_NUMBER(Pi)
2627
      P = Pr + (0.,1.)*Pi
2628

2629
      DO j = 1,s
2630
        DO k = 1,j-1
2631
          alpha(k) = dotprodfun(n, P(:,k), 1, P(:,j), 1 )
2632
          P(:,j) = P(:,j) - alpha(k)*P(:,k)
2633
        END DO
2634
        P(:,j) = P(:,j)/normfun(n,P(:,j),1)
2635
      END DO
2636
      kappa = 0.7d0
2637

2638
      ! Initialize local variables:
2639
      M = 0.0d0
2640
      om = 1.0d0
2641
      iter = 0
2642
      jj = 0
2643
      ii = 0
2644
      ! This concludes the initialisation phase
2645

2646
      ! Main iteration loop, build G-spaces:
2647

2648
      DO WHILE ( .NOT. Converged .AND. .NOT. Diverged )  ! start of iteration loop
2649

2650
        !!+++++++++++++++++++++++++++++++++++++++++++++++++++++++
2651
        ! Generate s vectors in G_j
2652
        !!+++++++++++++++++++++++++++++++++++++++++++++++++++++++
2653

2654
        ! New right-hand side for small system:
2655
        DO k = 1,s
2656
          f(k) = dotprodfun(n, P(:,k), 1, r, 1 )
2657
        END DO
2658

2659
        DO k = 1,s
2660

2661
          ! Update inner iteration counter
2662
          ii = ii + 1
2663

2664
          ! Compute new v
2665
          v = r 
2666
          IF ( jj > 0 ) THEN
2667

2668
            ! Solve small system (Note: M is lower triangular) and make v orthogonal to P:
2669
            DO i = k,s
2670
              gamma(i) = f(i)
2671
              DO j = k,i-1
2672
                gamma(i) = gamma(i) - M(i,j)*gamma(j)
2673
              END DO
2674
              gamma(i) = gamma(i)/M(i,i)
2675
              v = v - gamma(i)*G(:,i)
2676
            END DO
2677

2678
            ! Compute new U(:,k)
2679
            CALL pcondrsubr( t, v, ipar )
2680
            t = om*t
2681
            DO i = k,s
2682
              t = t + gamma(i)*U(:,i)
2683
            END DO
2684
            U(:,k) = t
2685

2686
          ELSE 
2687

2688
            ! Updates for the first s iterations (in G_0):
2689
            CALL pcondrsubr( U(:,k), v, ipar )
2690

2691
          END IF
2692

2693
          ! Compute new G(:,k), G(:,k) is in space G_j
2694
          CALL matvecsubr( U(:,k), G(:,k), ipar )
2695

2696
          ! Bi-Orthogonalise the new basis vectors: 
2697
          DO i = 1,s
2698
            mu(i) = dotprodfun(n, P(:,i), 1, G(:,k), 1 )
2699
          END DO
2700
          DO i = 1,k-1
2701
            alpha(i) = mu(i)
2702
            DO j = 1, i-1
2703
              alpha(i) = alpha(i) - M(i,j)*alpha(j)
2704
            END DO
2705
            alpha(i) = alpha(i)/M(i,i)
2706
            G(:,k) = G(:,k) - G(:,i)*alpha(i)
2707
            U(:,k) = U(:,k) - U(:,i)*alpha(i)
2708
            mu(k:s)  = mu(k:s)  - M(k:s,i)*alpha(i)
2709
          END DO
2710
          M(k:s,k) = mu(k:s)
2711

2712
          ! Break down?
2713
          IF ( ABS(M(k,k)) <= TINY(tol) ) THEN
2714
            Diverged = .TRUE.
2715
            EXIT
2716
          END IF
2717

2718
          ! Make r orthogonal to p_i, i = 1..k, update solution and residual 
2719
          beta(k) = f(k)/M(k,k)
2720
          r = r - beta(k)*G(:,k)
2721
          x = x + beta(k)*U(:,k)
2722

2723
          ! New f = P'*r (first k  components are zero)
2724
          IF ( k < s ) THEN
2725
            f(k+1:s)   = f(k+1:s) - beta(k)*M(k+1:s,k)
2726
          END IF
2727

2728
          ! Check for convergence
2729
          IF (UseStopCFun) THEN
2730
            errorind = stopcfun(x,b,r,ipar,dpar)
2731
          ELSE
2732
            normr = normfun(n,r,1)
2733
            errorind = normr/normb
2734
          END IF
2735
          iter = iter + 1
2736
          
2737
          IF( MOD(iter,OutputInterval) == 0) THEN
2738
            WRITE (*, '(I8, E11.4)') iter, errorind
2739
            CALL FLUSH(6)
2740
          END IF
2741

2742
          Converged = (errorind < Tol)
2743
          Diverged = (errorind > MaxTol) .OR. (errorind /= errorind)
2744
          IF ( Converged .OR. Diverged ) EXIT
2745
          IF (iter == MaxRounds) EXIT
2746

2747
        END DO ! Now we have computed s+1 vectors in G_j
2748
        IF ( Converged .OR. Diverged ) EXIT
2749
        IF (iter == MaxRounds) EXIT
2750

2751
        !!+++++++++++++++++++++++++++++++++++++++++++++++++++++++
2752
        ! Compute first residual in G_j+1
2753
        !!+++++++++++++++++++++++++++++++++++++++++++++++++++++++
2754

2755
        ! Update G-space counter
2756
        jj = jj + 1
2757

2758
        ! Compute first residual in G_j+1
2759
        ! Note: r is already perpendicular to P so v = r
2760

2761
        ! Preconditioning:
2762
        CALL pcondrsubr( v, r, ipar )
2763
        ! Matrix-vector multiplication:
2764
        CALL matvecsubr( v, t, ipar )
2765

2766
        ! Computation of a new omega
2767
        ! 'Maintaining the convergence':
2768
        nr = normfun(n,r,1)
2769
        nt = normfun(n,t,1)
2770
        tr = dotprodfun(n, t, 1, r, 1 )
2771
        rho = ABS(tr/(nt*nr))
2772
        om=tr/(nt*nt)
2773
        IF ( rho < kappa ) THEN
2774
          om = om*kappa/rho
2775
        END IF
2776

2777
        IF ( ABS(om) <= EPSILON(tol) ) THEN 
2778
          Diverged = .TRUE.
2779
          EXIT
2780
        END IF
2781

2782
        ! Update solution and residual
2783
        r = r - om*t 
2784
        x = x + om*v 
2785

2786
        ! Check for convergence
2787
        IF (UseStopCFun) THEN
2788
          errorind = stopcfun(x,b,r,ipar,dpar)
2789
        ELSE
2790
          normr = normfun(n,r,1)
2791
          errorind = normr/normb
2792
        END IF
2793
        iter = iter + 1
2794

2795
        IF( MOD(iter,OutputInterval) == 0) THEN
2796
          WRITE (*, '(I8, E11.4)') iter, errorind
2797
          CALL FLUSH(6)
2798
        END IF
2799

2800
        Converged = (errorind < Tol)
2801
        Diverged = (errorind > MaxTol) .OR. (errorind /= errorind)
2802
        IF (iter == MaxRounds) EXIT
2803

2804
      END DO ! end of while loop
2805

2806
    !----------------------------------------------------------
2807
    END SUBROUTINE ComplexIDRS
2808
    !----------------------------------------------------------
2809

2810
!--------------------------------------------------------------
2811
  END SUBROUTINE itermethod_z_idrs
2812
!--------------------------------------------------------------
2813

2814
END MODULE IterativeMethods
2815

2816
!> \} ElmerLib
2817

2818
Product

Resources

Company
