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), ALLOCATABLE :: P(:,:)
1602
      REAL(kind=dp), ALLOCATABLE :: G(:,:)
1603
      REAL(kind=dp), ALLOCATABLE :: U(:,:)
1604
      REAL(kind=dp), ALLOCATABLE :: r(:)
1605
      REAL(kind=dp), ALLOCATABLE :: v(:)
1606
      REAL(kind=dp), ALLOCATABLE :: t(:)
1607
      REAL(kind=dp), ALLOCATABLE :: M(:,:), f(:), mu(:)
1608
      REAL(kind=dp), ALLOCATABLE :: alpha(:), beta(:), gamma(:)
1609

1610
      REAL(kind=dp) :: om, tr, tr_s, tt
1611
      REAL(kind=dp) :: nr, nt, rho, kappa
1612

1613
      REAL(kind=dp), ALLOCATABLE :: r_s(:), x_s(:)
1614
      REAL(kind=dp) :: theta
1615
      
1616
      INTEGER :: iter                         ! number of iterations
1617
      INTEGER :: ii                           ! inner iterations index
1618
      INTEGER :: jj                           ! G-space index
1619
      REAL(kind=dp) :: normb, normr, errorind ! for tolerance check
1620
      INTEGER :: i,j,k,l                      ! loop counters
1621
!----------------------------------------------------------------------------------- 
1622

1623
      ALLOCATE( P(n,s), G(n,s), U(n,s), r(n), v(n), t(n), M(s,s), f(s), mu(s), alpha(s), beta(s), gamma(s))
1624
      
1625
      ! Compute initial residual
1626
      normb = normfun(n,b,1)
1627
      CALL C_matvec( x, t, ipar, matvecsubr )
1628
      r = b - t
1629
      
1630
      !-------------------------------------------------------------------
1631
      ! Check whether the initial guess satisfies the stopping criterion
1632
      !--------------------------------------------------------------------
1633
      IF (UseStopCFun) THEN
1634
        errorind = stopcfun(x,b,r,ipar,dpar)
1635
      ELSE
1636
        normr = normfun(n,r,1)
1637
        errorind = normr / normb
1638
      END IF
1639
      Converged = (errorind < Tol)
1640
      Diverged = (errorind > MaxTol) .OR. (errorind /= errorind)
1641

1642
      IF ( Converged .OR. Diverged ) RETURN
1643

1644
      U = 0.0d0
1645
      IF ( Smoothing ) THEN
1646
        ALLOCATE( r_s(n), x_s(n) )
1647
        x_s = x
1648
        r_s = r 
1649
      END IF
1650

1651
      ! Define P(n,s) and kappa
1652
#if 1
1653
      CALL RANDOM_NUMBER(P)
1654
#else
1655
      ! this is alternative generation of initial basis vectors
1656
      ! it is deterministic but not as good...
1657
      l = 0
1658
      k = 2        
1659
      DO j=1,s
1660
        DO i=1,n
1661
          P(i,j) = MODULO(i+l,k) / (1.0*(k-1)) 
1662
        END DO
1663
        l = k
1664
        k = 2*k + 1
1665
      END DO
1666
#endif
1667
              
1668
      DO j = 1,s
1669
        DO k = 1,j-1
1670
          alpha(k) = dotprodfun(n, P(:,k), 1, P(:,j), 1 )
1671
          P(:,j) = P(:,j) - alpha(k)*P(:,k)
1672
        END DO
1673
        P(:,j) = P(:,j)/normfun(n,P(:,j),1)
1674
      END DO
1675
      kappa = 0.7d0
1676

1677
      ! Initialize local variables:
1678
      M = 0.0d0
1679
      om = 1.0d0
1680
      iter = 0
1681
      jj = 0
1682
      ii = 0
1683

1684
      
1685
      ! This concludes the initialisation phase    
1686
      
1687
      ! Main iteration loop, build G-spaces:
1688
      DO WHILE ( (.NOT. Converged) .AND. (.NOT. Diverged) ) 
1689

1690
        !!+++++++++++++++++++++++++++++++++++++++++++++++++++++++
1691
        ! Generate s vectors in G_j
1692
        !!+++++++++++++++++++++++++++++++++++++++++++++++++++++++
1693

1694
        ! New right-hand side for small system:
1695
        DO k = 1,s
1696
          f(k) = dotprodfun(n, P(:,k), 1, r, 1 )
1697
        END DO
1698

1699
        DO k = 1,s
1700

1701
          ! Update inner iteration counter
1702
          ii = ii + 1
1703

1704
          ! Compute new v
1705
          v = r
1706
          IF ( jj > 0 ) THEN
1707

1708
            ! Solve small system (Note: M is lower triangular) and make v orthogonal to P:
1709
            DO i = k,s
1710
              gamma(i) = f(i)
1711
              DO j = k,i-1
1712
                gamma(i) = gamma(i) - M(i,j)*gamma(j)
1713
              END DO
1714
              gamma(i) = gamma(i)/M(i,i)
1715
              v = v - gamma(i)*G(:,i)
1716
            END DO
1717

1718
            ! Compute new U(:,k)
1719
            CALL C_rpcond( t, v, ipar, pcondrsubr ) 
1720
            t = om*t
1721
            DO i = k,s
1722
              t = t + gamma(i)*U(:,i)
1723
            END DO
1724
            U(:,k) = t
1725

1726
          ELSE
1727

1728
            ! Updates for the first s iterations (in G_0):
1729
            CALL C_rpcond( U(:,k), v, ipar, pcondrsubr )
1730

1731
          END IF
1732

1733
          ! Compute new G(:,k), G(:,k) is in space G_j
1734
          CALL C_matvec( U(:,k), G(:,k), ipar, matvecsubr )
1735

1736
          ! Bi-Orthogonalise the new basis vectors:
1737
          DO i = 1,s
1738
            mu(i) = dotprodfun(n, P(:,i), 1, G(:,k), 1 )
1739
          END DO
1740
          DO i = 1,k-1
1741
            alpha(i) = mu(i)
1742
            DO j = 1, i-1
1743
              alpha(i) = alpha(i) - M(i,j)*alpha(j)
1744
            END DO
1745
            alpha(i) = alpha(i)/M(i,i)
1746
            G(:,k) = G(:,k) - G(:,i)*alpha(i)
1747
            U(:,k) = U(:,k) - U(:,i)*alpha(i)
1748
            mu(k:s)  = mu(k:s)  - M(k:s,i)*alpha(i)
1749
          END DO
1750
          M(k:s,k) = mu(k:s)
1751

1752
          ! Break down?
1753
          IF ( ABS(M(k,k)) <= TINY(tol) ) THEN
1754
            Diverged = .TRUE.
1755
            EXIT
1756
          END IF
1757

1758
          ! Make r orthogonal to p_i, i = 1..k, update solution and residual
1759
          beta(k) = f(k)/M(k,k)
1760
          r = r - beta(k)*G(:,k)
1761
          x = x + beta(k)*U(:,k)
1762
          
1763
          ! New f = P'*r (first k  components are zero)
1764
          IF ( k < s ) THEN
1765
            f(k+1:s)   = f(k+1:s) - beta(k)*M(k+1:s,k)
1766
          END IF
1767

1768
          IF (Smoothing) THEN
1769
            t = r_s - r
1770
            tr_s = dotprodfun(n, t, 1, r_s, 1 )
1771
            tt = dotprodfun(n, t, 1, t, 1 )
1772
            theta = tr_s / tt
1773
            
1774
            r_s = r_s - theta * t
1775
            x_s = x_s - theta * (x_s - x)
1776
          END IF
1777

1778
          ! Check for convergence
1779
          iter = iter + 1
1780
          IF (UseStopCFun) THEN
1781
            IF (Smoothing) THEN
1782
              errorind = stopcfun(x_s,b,r_s,ipar,dpar)
1783
            ELSE
1784
              errorind = stopcfun(x,b,r,ipar,dpar)
1785
            END IF
1786
          ELSE
1787
            IF (Smoothing) THEN
1788
              normr = normfun(n,r_s,1)  
1789
            ELSE
1790
              normr = normfun(n,r,1)
1791
            END IF
1792
            errorind = normr/normb
1793
          END IF
1794

1795
          IF( MOD(iter,OutputInterval) == 0) THEN
1796
            WRITE (*, '(I8, E11.4)') iter, errorind
1797
          END IF
1798

1799
          Converged = (errorind < Tol)
1800
          Diverged = (errorind > MaxTol) .OR. (errorind /= errorind)
1801
          IF ( Converged .OR. Diverged ) EXIT
1802
          IF (iter == MaxRounds) EXIT
1803
          
1804
         
1805
        END DO ! Now we have computed s+1 vectors in G_j
1806
        IF ( Converged .OR. Diverged ) EXIT
1807
        IF (iter == MaxRounds) EXIT
1808

1809
        !!+++++++++++++++++++++++++++++++++++++++++++++++++++++++
1810
        ! Compute first residual in G_j+1
1811
        !!+++++++++++++++++++++++++++++++++++++++++++++++++++++++
1812

1813
        ! Update G-space counter
1814
        jj = jj + 1
1815

1816
        ! Compute first residual in G_j+1
1817
        ! Note: r is already perpendicular to P so v = r
1818

1819
        ! Preconditioning:
1820
        CALL C_rpcond( v, r, ipar, pcondrsubr )
1821
        ! Matrix-vector multiplication:
1822
        CALL C_matvec( v, t, ipar, matvecsubr )
1823

1824
        ! Computation of a new omega
1825
        ! 'Maintaining the convergence':
1826
        nr = normfun(n,r,1)
1827
        nt = normfun(n,t,1)
1828
        tr = dotprodfun(n, t, 1, r, 1 )
1829
        rho = ABS(tr/(nt*nr))
1830
        om=tr/(nt*nt)
1831
        IF ( rho < kappa ) THEN
1832
          om = om*kappa/rho
1833
        END IF
1834

1835
        IF ( ABS(om) <= EPSILON(tol) ) THEN
1836
          Diverged = .TRUE.
1837
          EXIT
1838
        END IF
1839

1840
        ! Update solution and residual
1841
        r = r - om*t
1842
        x = x + om*v
1843
        
1844
        IF (Smoothing) THEN
1845
          t = r_s - r
1846
          tr_s = dotprodfun(n, t, 1, r_s, 1 )
1847
          tt = dotprodfun(n, t, 1, t, 1 )
1848
          theta = tr_s / tt
1849
          r_s = r_s - theta * t
1850
          x_s = x_s - theta * (x_s - x)
1851
        END IF
1852

1853
        ! Check for convergence
1854
        iter = iter + 1
1855
        IF (UseStopCFun) THEN
1856
          IF (Smoothing) THEN
1857
            errorind = stopcfun(x_s,b,r_s,ipar,dpar)
1858
          ELSE
1859
            errorind = stopcfun(x,b,r,ipar,dpar)
1860
          END IF
1861
        ELSE
1862
          IF (Smoothing) THEN
1863
            normr = normfun(n,r_s,1)  
1864
          ELSE
1865
            normr = normfun(n,r,1)
1866
          END IF
1867
          errorind = normr/normb
1868
        END IF
1869
        
1870
        IF( MOD(iter,OutputInterval) == 0) THEN
1871
          WRITE (*, '(I8, E11.4)') iter, errorind
1872
          CALL FLUSH(6)
1873
        END IF
1874

1875
        IF( Robust ) THEN
1876
          ! Always store the best solution so far (with some small margin)
1877
          IF( errorInd < RobustStep * BestNorm ) THEN
1878
            BestIter = iter
1879
            BestNorm = errorInd
1880
            IF (Smoothing) THEN
1881
              Bestx = x_s
1882
            ELSE
1883
              Bestx = x
1884
            END IF
1885
            BadIterCount = 0
1886
          ELSE
1887
            BadIterCount = BadIterCount + 1
1888
          END IF
1889

1890
          ! If we have diverged too much and have found already a good candidate, then take it
1891
          IF( BestNorm <  RobustTol .AND. &
1892
              ( errorInd > RobustMaxTol .OR. BadIterCount > MaxBadIter ) ) THEN
1893
            EXIT
1894
          END IF
1895
          
1896
        END IF
1897
                        
1898
        Converged = (errorind < Tol)
1899
        Diverged = (errorind > MaxTol) .OR. (errorind /= errorind)
1900
        IF (iter == MaxRounds) EXIT
1901
      END DO ! end of while loop
1902

1903
      IF( Smoothing ) x = x_s
1904
      
1905
      IF( Robust ) THEN
1906
        IF( BestNorm < RobustTol ) THEN
1907
          Converged = .TRUE.
1908
        END IF
1909
        IF( BestNorm < errorInd ) THEN
1910
          x = Bestx
1911
        END IF        
1912
        IF(OutputInterval /= HUGE(OutputInterval)) THEN
1913
          WRITE(*,'(A,I8,E11.4,I8,E11.4)') 'Idrs robust: ',&
1914
              iter, BestNorm, BestIter, errorind
1915
          CALL FLUSH(6)
1916
        END IF
1917
      ELSE
1918
        IF(OutputInterval /= HUGE(OutputInterval)) THEN
1919
          WRITE(*,'(A,I8,E11.4)') 'Idrs: ',&
1920
              iter, errorind
1921
          CALL FLUSH(6)
1922
        END IF
1923
      END IF
1924
      
1925
    !----------------------------------------------------------
1926
    END SUBROUTINE RealIDRS
1927
    !----------------------------------------------------------
1928

1929
!--------------------------------------------------------------
1930
  END SUBROUTINE itermethod_idrs
1931
!--------------------------------------------------------------
1932

1933
!---------------------------------------------------------------------------------
1934
!>  This routine solves real linear system Ax = b with inequality constraint using
1935
!>  the MPRGP algorithm (Modiﬁed Proportioning and Reduced Gradient Projection).
1936
!---------------------------------------------------------------------------------
1937
  SUBROUTINE itermethod_mprgp( xvec, rhsvec, &
1938
      ipar, dpar, work, matvecsubr, pcondlsubr, &
1939
      pcondrsubr, dotprodfun, normfun, stopcfun )
1940
!------------------------------------------------------------------------------
1941
    USE huti_interfaces
1942
    IMPLICIT NONE
1943
    PROCEDURE( mv_iface_d ), POINTER :: matvecsubr
1944
    PROCEDURE( pc_iface_d ), POINTER :: pcondlsubr
1945
    PROCEDURE( pc_iface_d ), POINTER :: pcondrsubr
1946
    PROCEDURE( dotp_iface_d ), POINTER :: dotprodfun
1947
    PROCEDURE( norm_iface_d ), POINTER :: normfun
1948
    PROCEDURE( stopc_iface_d ), POINTER :: stopcfun
1949

1950
    INTEGER, DIMENSION(HUTI_IPAR_DFLTSIZE) :: ipar
1951
    REAL(KIND=dp), TARGET, DIMENSION(HUTI_NDIM) :: xvec, rhsvec
1952
    ! DOUBLE PRECISION, DIMENSION(HUTI_NDIM), TARGET :: xvec, rhsvec
1953
    DOUBLE PRECISION, DIMENSION(HUTI_DPAR_DFLTSIZE) :: dpar
1954
    REAL(KIND=dp), DIMENSION(HUTI_WRKDIM,HUTI_NDIM) :: work
1955
    ! DOUBLE PRECISION, DIMENSION(HUTI_WRKDIM,HUTI_NDIM) :: work
1956
    INTEGER :: ndim,i,j,k
1957
    INTEGER :: Rounds, OutputInterval, s
1958
    REAL(KIND=dp) :: MinTol, MaxTol, Residual
1959
    LOGICAL :: Converged, Diverged, UseStopCFun
1960
    INTEGER :: ncg, ne, np, iters
1961
    REAL(KIND=dp) :: final_norm_gp
1962
    
1963
    INTEGER :: n
1964

1965
    ! MPRGP parameters (passed from calling routine)
1966
    REAL(KIND=dp) :: Gamma, TolFactor
1967
    LOGICAL :: adapt
1968
    INTEGER :: bound
1969

1970
    TYPE(Matrix_t), POINTER :: A, Aser
1971

1972
    REAL(KIND=dp), POINTER :: x(:),b(:),c(:),cser(:)
1973

1974
    A => GlobalMatrix    
1975
    ndim = HUTI_NDIM
1976

1977
    x => xvec
1978
    b => rhsvec
1979

1980
    Rounds = HUTI_MAXIT
1981
    MinTol = HUTI_TOLERANCE
1982
    MaxTol = HUTI_MAXTOLERANCE
1983
    OutputInterval = HUTI_DBUGLVL 
1984
    s = HUTI_IDRS_S
1985
    UseStopCFun = HUTI_STOPC == HUTI_USUPPLIED_STOPC
1986
    
1987
    Gamma = HUTI_MPRGP_GAMMA
1988
    Adapt = ( HUTI_MPRGP_ADAPT > 0 )
1989
    TolFactor = HUTI_MPRGP_TOLFACTOR
1990

1991
    IF (ASSOCIATED(CurrentModel % Solver % Variable % LowerLimit)) THEN
1992
      bound = 0
1993
      cser => CurrentModel % Solver % Variable % LowerLimit
1994
    ELSE IF (ASSOCIATED(CurrentModel % Solver % Variable % UpperLimit)) THEN
1995
      bound = 1
1996
      cser => CurrentModel % Solver % Variable % UpperLimit
1997
    ELSE
1998
      CALL Fatal('itermethod_mprgp','No limits found')
1999
    END IF
2000

2001
    ! Scaling here is a little dirty, but this is still under development...
2002
    ! The serial matrix and limiter have different size than the parallel ones.
2003
    Aser => CurrentModel % Solver % Matrix
2004
    IF( ASSOCIATED(Aser % DiagScaling ) ) THEN
2005
      cser(1:ndim) = cser(1:ndim) / (Aser % DiagScaling * Aser % RhsScaling) 
2006
    END IF
2007
    IF( ParEnv % PEs == 1) THEN
2008
      c => cser
2009
    ELSE
2010
      ALLOCATE(c(ndim))
2011
      j = 0
2012
      ! We need to map the bigger limiter to the parallel one where only owned dofs are considered.
2013
      DO i=1,Aser % NumberOfRows
2014
        IF (Aser % ParallelInfo % NeighbourList(i) % Neighbours(1) == ParEnv % Mype ) THEN
2015
          j = j + 1
2016
          c(j) = cser(i)
2017
        END IF
2018
      END DO
2019
    END IF
2020
      
2021
    
2022
    Converged = .FALSE.
2023
    Diverged = .FALSE.
2024
    n = ndim
2025
    CALL MPRGP(ndim, x, b, c, MinTol, Rounds, Gamma, adapt, bound, TolFactor, &
2026
        ncg, ne, np, iters, Converged, final_norm_gp )
2027

2028
    CALL Info('MPRGPSolver','MPRGP finished: iters='//I2S(iters)// &
2029
          ', ncg='//I2S(ncg)//', ne='//I2S(ne)//', np='//I2S(np), Level=10)
2030

2031
    IF(Converged) HUTI_INFO = HUTI_CONVERGENCE
2032
    IF(Diverged) HUTI_INFO = HUTI_DIVERGENCE
2033
    IF ( (.NOT. Converged) .AND. (.NOT. Diverged) ) HUTI_INFO = HUTI_MAXITER
2034

2035

2036
    IF( ASSOCIATED(Aser % DiagScaling ) ) THEN
2037
      cser(1:ndim) = cser(1:ndim) * (Aser % DiagScaling * Aser % RhsScaling) 
2038
    END IF
2039
    IF( ParEnv % PEs > 1) DEALLOCATE(c)
2040

2041
    
2042
  CONTAINS
2043

2044
!-----------------------------------------------------------------------------------
2045
!   Implementation of the MPRGP algorithm. 
2046
!   The subroutine MPRGP was written by D. Reeves during an internship at CSC.
2047
!----------------------------------------------------------------------------------- 
2048
    SUBROUTINE MPRGP(n, x, b, c, epsr, maxit, Gamma, adapt, bound, TolFactor, &
2049
                      ncg, ne, np, iters, converged, final_norm_gp)
2050
!----------------------------------------------------------------------------------- 
2051
      ! ---------------------------
2052
      ! Arguments
2053
      ! ---------------------------
2054
      INTEGER, INTENT(IN) :: n
2055
      REAL(KIND=dp), INTENT(INOUT) :: x(n)    ! initial guess in, final solution out
2056
      REAL(KIND=dp), INTENT(IN) :: b(n), c(n) ! rhs and bound
2057
      REAL(KIND=dp), INTENT(IN) :: epsr
2058
      INTEGER, INTENT(IN) :: maxit
2059
      REAL(KIND=dp), INTENT(IN) :: Gamma
2060
      LOGICAL, INTENT(IN) :: adapt
2061
      INTEGER, INTENT(IN) :: bound      ! 'lower' or 'upper'
2062
      REAL(KIND=dp), INTENT(IN) :: TolFactor
2063
      INTEGER, INTENT(OUT) :: ncg, ne, np, iters
2064
      LOGICAL, INTENT(OUT) :: converged
2065
      REAL(KIND=dp), INTENT(OUT) :: final_norm_gp
2066
      INTEGER :: ierr, comm, n_beyond
2067
      
2068
      INTEGER :: itl
2069
      REAL(KIND=dp) :: normv
2070
      
2071
      ! ---------------------------
2072
      ! Local declarations (all here)
2073
      ! ---------------------------
2074
      REAL(KIND=dp), ALLOCATABLE :: g(:), gf(:), gc(:), gr(:), gp(:)
2075
      REAL(KIND=dp), ALLOCATABLE :: z(:), p(:), Ap(:), yy(:), D(:), Agr(:)
2076
      LOGICAL, ALLOCATABLE :: J(:)
2077
      LOGICAL, ALLOCATABLE :: p_mask(:)
2078
      INTEGER :: bs
2079
      REAL(KIND=dp) :: lAl, alpha, a_f
2080
      INTEGER :: i, allocstat
2081
      REAL(KIND=dp) :: rtp, pAp, acg, beta
2082
      REAL(KIND=dp) :: grg, grAgr
2083
      REAL(KIND=dp) :: eps_local
2084
      TYPE(Matrix_t), POINTER :: MatA
2085
      REAL(KIND=dp) :: tol
2086
      
2087
      REAL(KIND=dp), ALLOCATABLE :: v(:),w(:)
2088

2089
      comm = ParEnv % ActiveComm
2090

2091
      ! ---------------------------
2092
      ! Allocate scratch vectors
2093
      ! ---------------------------
2094
      ALLOCATE(g(n), gf(n), gc(n), gr(n), gp(n), z(n), p(n), Ap(n), yy(n), J(n), Agr(n), STAT=allocstat)
2095
      IF (allocstat /= 0) THEN
2096
        CALL Fatal('itermethod_mprgp','Allocation failed for scratch vectors')
2097
      END IF
2098

2099
      ! ---------------------------
2100
      ! Preliminaries
2101
      ! ---------------------------
2102
      eps_local = EPSILON(1.0_dp)
2103

2104
      ! set bound sign
2105
      IF (bound == 1) THEN
2106
        bs = -1
2107
      ELSE
2108
        bs = 1
2109
      END IF
2110

2111

2112
      ! ---------------------------
2113
      ! Initialization (g = A*x - b)
2114
      ! ---------------------------
2115
      CALL matvecsubr( x, g, ipar )
2116

2117
      g = g - b
2118

2119
      tol = 1.0e-12_dp
2120
      IF (bs == 1) THEN
2121
        J = (x > c + tol) ! J is free set
2122
      ELSE
2123
        J = (x < c - tol)
2124
      END IF
2125
      gf = MERGE(g, 0.0_dp, J)
2126

2127
      IF (bs == 1) THEN
2128
        WHERE (.NOT. J)
2129
          gc = MIN(g, 0.0_dp)
2130
        ELSEWHERE
2131
          gc = 0.0_dp
2132
        END WHERE
2133
      ELSE
2134
        WHERE (.NOT. J)  ! WHERE COMMAND FOR vECTORS
2135
          gc = MAX(g, 0.0_dp)
2136
        ELSEWHERE
2137
          gc = 0.0_dp
2138
        END WHERE
2139
      END IF
2140

2141
      ! estimate matrix norm ||A|| with a small power iteration
2142
      CALL estimate_matrix_norm(n, 10, lAl)
2143
!      lal = 1.0_dp
2144
      IF (lAl <= 0.0_dp) lAl = 1.0_dp
2145
      alpha = 1.0_dp / lAl
2146

2147
      ! reduced free gradient gr
2148
      IF (bs == 1) THEN
2149
        gr = MERGE(MIN(lAl * (x - c), gf), 0.0_dp, J)
2150
      ELSE
2151
        gr = MERGE(MAX(lAl * (x - c), gf), 0.0_dp, J)
2152
      END IF
2153

2154
      gp = gf + gc
2155
      ! preconditioning: z = M^{-1} * g on free set
2156
      z = g
2157
      CALL pcondrsubr( g, z, ipar )
2158
      WHERE (.NOT. J)
2159
        z = 0.0_dp
2160
      END WHERE
2161
      p = z
2162

2163
      ! counters
2164
      ncg = 0
2165
      ne = 0
2166
      np = 0
2167
      iters = 0
2168
      converged = .FALSE.
2169

2170
      ! ---------------------------
2171
      ! Main loop
2172
      ! ---------------------------
2173
      DO WHILE ( normfun(n, gp, 1) > epsr .AND. iters < maxit )
2174

2175
        iters = iters + 1
2176
        IF ( dotprodfun(n, gc, 1, gc, 1) <= (Gamma**2) * dotprodfun(n, gr, 1, gf, 1) ) THEN
2177
          ! CG step
2178
          CALL matvecsubr( p, Ap, ipar )        
2179
          rtp = dotprodfun(n, z, 1, g, 1) ! residual * p
2180
          pAp = dotprodfun(n, p, 1, Ap, 1)
2181

2182
          IF (ABS(pAp) < eps_local) THEN
2183
            ! We don't want to divide by zero, if norm of gp is small enough, we can stop
2184
            IF (normfun(n, gp, 1) < TolFactor * epsr) THEN
2185
              converged = .TRUE.
2186
              CALL Info('itermethod_mprgp','MPRGP converged')
2187
              EXIT
2188
            END IF
2189
            CALL Fatal('itermethod_mprgp','p''*A*p nearly zero, stopping')
2190
          END IF
2191

2192
          acg = rtp / pAp
2193
          yy = x - acg * p
2194

2195
          n_beyond = COUNT(bs * yy(:) < bs * c(:))
2196
          IF(ParEnv % PEs > 1) THEN
2197
            CALL MPI_ALLREDUCE( MPI_IN_PLACE, n_beyond, 1, MPI_INTEGER, MPI_SUM, comm, ierr )              
2198
          END IF
2199

2200
          IF (n_beyond == 0) THEN
2201
            ! accept full CG step
2202
            x = yy
2203
            g = g - acg * Ap
2204
            IF (bs == 1) THEN
2205
              J = (x > c + tol)
2206
            ELSE
2207
              J = (x < c - tol)
2208
            END IF
2209

2210
            ! precondition
2211
            z = g
2212
            CALL pcondrsubr( g, z, ipar )
2213
            WHERE (.NOT. J)
2214
              z = 0.0_dp
2215
            END WHERE
2216

2217
            beta = dotprodfun(n, z, 1, Ap, 1) / pAp
2218
            p = z - beta * p
2219

2220
            ! update gf,gc,gr,gp
2221
            gf = MERGE(g, 0.0_dp, J)
2222

2223
            IF (bs == 1) THEN
2224
              gc = 0.0_dp
2225
              WHERE (.NOT. J)
2226
                gc = MIN(g, 0.0_dp)
2227
              END WHERE
2228
              gr = MERGE(MIN(lAl * (x - c), gf), 0.0_dp, J)
2229
            ELSE
2230
              gc = 0.0_dp
2231
              WHERE (.NOT. J)
2232
                gc = MAX(g, 0.0_dp)
2233
              END WHERE
2234
              gr = MERGE(MAX(lAl * (x - c), gf), 0.0_dp, J)
2235
            END IF
2236
            gp = gf + gc
2237
            ncg = ncg + 1
2238
          ELSE
2239
            ! expansion step (feasible step)
2240
            p_mask = (bs * p > 0.0_dp) .AND. J ! indexes where p is moving towards the bound
2241
            a_f = MINVAL((x-c) / p,p_mask)
2242
            
2243
            IF(ParEnv % PEs > 1) THEN
2244
              CALL MPI_ALLREDUCE( MPI_IN_PLACE, a_f, 1, MPI_DOUBLE_PRECISION, MPI_MIN, comm, ierr )              
2245
            END IF
2246
              
2247
            IF (a_f < 0.0_dp) a_f = 0.0_dp
2248
            ! halfstep
2249
            IF (bs == 1) THEN
2250
                x = MAX(x - a_f * p, c)
2251
            ELSE
2252
                x = MIN(x - a_f * p, c)
2253
            END IF
2254

2255
            IF (bs == 1) THEN
2256
              J = (x > c + tol)
2257
            ELSE
2258
              J = (x < c - tol)
2259
            END IF
2260
            g = g - a_f * Ap
2261

2262
            ! adaptive alpha
2263
            IF (adapt) THEN
2264
              CALL matvecsubr( gr, Agr, ipar )        
2265
              !CALL C_matvec(gr, Agr, ipar, matvecsubr)
2266
              grg = dotprodfun(n, gr, 1, g, 1)
2267
              grAgr = dotprodfun(n, gr, 1, Agr, 1)
2268
              IF (ABS(grAgr) < eps_local) THEN
2269
                alpha = 1.0_dp / lAl
2270
              ELSE
2271
                alpha = grg / grAgr
2272
                IF (alpha <= 0.0_dp .OR. alpha > 1.0_dp / lAl) alpha = 1.0_dp / lAl
2273
              END IF
2274
            END IF
2275

2276
            IF (bs == 1) THEN
2277
              WHERE (J)
2278
                x = MAX(x - alpha * g, c)
2279
              END WHERE
2280
            ELSE
2281
              WHERE (J)
2282
                x = MIN(x - alpha * g, c)
2283
              END WHERE
2284
            END IF
2285

2286
            CALL matvecsubr( x, g, ipar )
2287
            g = g - b
2288

2289
            ! recompute z and p
2290
            z = g
2291
            CALL pcondrsubr( g, z, ipar )
2292
            WHERE (.NOT. J)
2293
              z = 0.0_dp
2294
            END WHERE
2295
            p = z
2296

2297
            ! update gf,gc,gr,gp
2298
            gf = MERGE(g, 0.0_dp, J)
2299

2300
            IF (bs == 1) THEN
2301
              gc = 0.0_dp
2302
              WHERE (.NOT. J)
2303
                gc = MIN(g, 0.0_dp)
2304
              END WHERE
2305
              gr = MERGE(MIN(lAl * (x - c), gf), 0.0_dp, J)
2306
            ELSE
2307
              gc = 0.0_dp
2308
              WHERE (.NOT. J)
2309
                gc = MAX(g, 0.0_dp)
2310
              END WHERE
2311
              gr = MERGE(MAX(lAl * (x - c), gf), 0.0_dp, J)
2312
            END IF
2313

2314
            gp = gf + gc
2315
            ne = ne + 1
2316
          END IF
2317

2318
        ELSE
2319
          ! proportioning step
2320
          CALL matvecsubr( gc, Ap, ipar )    ! A*gc
2321
          pAp = dotprodfun(n, gc, 1, Ap, 1)
2322
          IF (ABS(pAp) < eps_local) THEN
2323
            IF (normfun(n, gp, 1) < TolFactor * epsr) THEN
2324
              converged = .TRUE.
2325
              CALL Info('itermethod_mprgp','MPRGP converged')
2326
              EXIT
2327
            END IF
2328

2329
            CALL FATAL('itermethod_mprgp','denominator gc*A*gc in proportioning step nearly zero, stopping')
2330
            EXIT
2331
          END IF
2332
          acg = dotprodfun(n, gc, 1, g, 1) / pAp
2333

2334
          x = x - acg * gc
2335

2336
          ! safeguard if we end up outside the bound
2337
          IF (bs == 1) THEN
2338
            x = MAX(x, c)
2339
          ELSE
2340
            x = MIN(x, c)
2341
          END IF
2342

2343
          IF (bs == 1) THEN
2344
            J = (x > c + tol)
2345
          ELSE
2346
            J = (x < c - tol)
2347
          END IF
2348
          g = g - acg * Ap
2349
          
2350
          z = g
2351
          CALL pcondrsubr( g, z, ipar )
2352
          WHERE (.NOT. J)
2353
            z = 0.0_dp
2354
          END WHERE
2355
          p = z
2356

2357
          ! update gf,gc,gr,gp
2358
          gf = MERGE(g, 0.0_dp, J)
2359

2360
          IF (bs == 1) THEN
2361
            gc = 0.0_dp
2362
            WHERE (.NOT. J)
2363
              gc = MIN(g, 0.0_dp)
2364
            END WHERE
2365
            gr = MERGE(MIN(lAl * (x - c), gf), 0.0_dp, J)
2366
          ELSE
2367
            gc = 0.0_dp
2368
            WHERE (.NOT. J)
2369
              gc = MAX(g, 0.0_dp)
2370
            END WHERE
2371
            gr = MERGE(MAX(lAl * (x - c), gf), 0.0_dp, J)
2372
          END IF
2373

2374
          gp = gf + gc
2375
          np = np + 1
2376
        END IF
2377

2378
      END DO  ! main loop
2379

2380
      final_norm_gp = normfun(n, gp, 1)
2381
      
2382
      IF (.NOT. converged) THEN ! can be set as converged in the main loop
2383
        converged = (final_norm_gp <= epsr)
2384
      END IF
2385

2386
      ! cleanup
2387
      IF (ALLOCATED(D)) DEALLOCATE(D)
2388
      DEALLOCATE(g, gf, gc, gr, gp, z, p, Ap, yy, J)
2389

2390
    !----------------------------------------------------------
2391
    END SUBROUTINE MPRGP
2392
    !----------------------------------------------------------
2393

2394
    ! Estimate matrix norm using power method
2395
    SUBROUTINE estimate_matrix_norm(nloc, niter, out_norm)
2396
      INTEGER, INTENT(IN) :: nloc, niter
2397
      REAL(KIND=dp), INTENT(OUT) :: out_norm
2398
      TYPE(Matrix_t), POINTER :: A
2399
      REAL(KIND=dp), ALLOCATABLE :: v(:), w(:)
2400
      INTEGER :: itl, allocstat
2401
      REAL(KIND=dp) :: normv
2402

2403
      ALLOCATE(v(nloc), w(nloc), STAT=allocstat)
2404
      IF (allocstat /= 0) THEN
2405
        CALL Fatal('estimate_matrix_norm','alloc fail')
2406
      END IF
2407

2408
      v = 1.0_dp / REAL(nloc, dp)
2409

2410
      DO itl = 1, niter
2411
        !CALL C_matvec(v, w, ipar, matvecsubr)
2412
        CALL matvecsubr(v, w, ipar)
2413
        normv = normfun(nloc, w, 1)
2414
        IF (normv <= 0.0_dp) EXIT
2415
        v = w / normv
2416
      END DO
2417

2418
      !CALL matvecsubr(v, w, ipar)
2419
      CALL matvecsubr(v, w, ipar)
2420
      out_norm = normfun(nloc, w, 1)
2421
      DEALLOCATE(v, w)
2422
    END SUBROUTINE estimate_matrix_norm
2423

2424
!--------------------------------------------------------------
2425
  END SUBROUTINE itermethod_mprgp
2426
!--------------------------------------------------------------
2427

2428
!------------------------------------------------------------------------------
2429
!> This routine provides the complex version to the GCR linear solver.
2430
!------------------------------------------------------------------------------
2431
 SUBROUTINE itermethod_z_gcr( xvec, rhsvec, &
2432
      ipar, dpar, work, matvecsubr, pcondlsubr, &
2433
      pcondrsubr, dotprodfun, normfun, stopcfun )
2434
!------------------------------------------------------------------------------
2435
    USE huti_interfaces
2436
    IMPLICIT NONE
2437
    PROCEDURE( mv_iface_z ), POINTER :: matvecsubr
2438
    PROCEDURE( pc_iface_z ), POINTER :: pcondlsubr
2439
    PROCEDURE( pc_iface_z ), POINTER :: pcondrsubr
2440
    PROCEDURE( dotp_iface_z ), POINTER :: dotprodfun
2441
    PROCEDURE( norm_iface_z ), POINTER :: normfun
2442
    PROCEDURE( stopc_iface_z ), POINTER :: stopcfun
2443

2444
    INTEGER, DIMENSION(HUTI_IPAR_DFLTSIZE) :: ipar
2445
    COMPLEX(KIND=dp), TARGET, DIMENSION(HUTI_NDIM) :: xvec, rhsvec
2446
    DOUBLE PRECISION, DIMENSION(HUTI_DPAR_DFLTSIZE) :: dpar
2447
    COMPLEX(KIND=dp), DIMENSION(HUTI_WRKDIM,HUTI_NDIM) :: work
2448

2449
    COMPLEX(KIND=dp) :: y(HUTI_NDIM),f(HUTI_NDIM)
2450
    INTEGER :: ndim, RestartN, i
2451
    INTEGER :: Rounds, OutputInterval
2452
    REAL(KIND=dp) :: MinTol, MaxTol, Residual
2453
    LOGICAL :: Converged, Diverged, UseStopCFun
2454

2455
    CALL Info('Itermethod_z_gcr','Starting GCR iteration',Level=25)
2456
    
2457
    ndim = HUTI_NDIM
2458
    Rounds = HUTI_MAXIT
2459
    MinTol = HUTI_TOLERANCE
2460
    MaxTol = HUTI_MAXTOLERANCE
2461
    OutputInterval = HUTI_DBUGLVL
2462
    RestartN = HUTI_GCR_RESTART 
2463
    UseStopCFun = HUTI_STOPC == HUTI_USUPPLIED_STOPC
2464
    
2465
    Converged = .FALSE.
2466
    Diverged = .FALSE.
2467
    
2468
    !----------------------------------------------------------------------------
2469
    ! Transform the solution vector and the right-hand side vector to 
2470
    ! complex-valued vectors y and f
2471
    !---------------------------------------------------------------------------
2472
    DO i=1,ndim
2473
      y(i)=xvec(i)
2474
      f(i)=rhsvec(i)
2475
    END DO
2476
       
2477
    CALL GCR_Z(ndim, GlobalMatrix, y, f, Rounds, MinTol, MaxTol, Residual, &
2478
        Converged, Diverged, OutputInterval, RestartN )
2479

2480
    IF(Converged) HUTI_INFO = HUTI_CONVERGENCE
2481
    IF(Diverged) HUTI_INFO = HUTI_DIVERGENCE
2482
    IF ( (.NOT. Converged) .AND. (.NOT. Diverged) ) HUTI_INFO = HUTI_MAXITER
2483
   
2484
    DO i=1,ndim
2485
      xvec(i) = y(i)
2486
    END DO
2487

2488
  CONTAINS 
2489
    
2490
    
2491
!------------------------------------------------------------------------------  
2492
    SUBROUTINE GCR_Z( n, A, x, b, Rounds, MinTolerance, MaxTolerance, Residual, &
2493
        Converged, Diverged, OutputInterval, m) 
2494
!------------------------------------------------------------------------------
2495
      TYPE(Matrix_t), POINTER :: A
2496
      INTEGER :: Rounds
2497
      COMPLEX(KIND=dp) :: x(n),b(n)
2498
      LOGICAL :: Converged, Diverged
2499
      REAL(KIND=dp) :: MinTolerance, MaxTolerance, Residual
2500
      INTEGER :: n, OutputInterval, m
2501
!------------------------------------------------------------------------------      
2502
      REAL(KIND=dp) :: bnorm,rnorm
2503
      COMPLEX(KIND=dp), ALLOCATABLE :: R(:)
2504
      COMPLEX(KIND=dp), ALLOCATABLE :: S(:,:), V(:,:), T1(:), T2(:)
2505
      INTEGER :: i,j,k,allocstat
2506
      COMPLEX(KIND=dp) :: beta, czero
2507
      REAL(KIND=dp) :: alpha, trueresnorm, normerr
2508
      COMPLEX(KIND=dp), ALLOCATABLE :: trueres(:)
2509
!------------------------------------------------------------------------------
2510
            
2511
      ALLOCATE( R(n), T1(n), T2(n), trueres(n), STAT=allocstat )
2512
      IF( allocstat /= 0) &
2513
          CALL Fatal('GCR_Z','Failed to allocate memory of size: '//I2S(n))
2514
      IF ( m > 1 ) THEN        
2515
         ALLOCATE( S(n,m-1), V(n,m-1), STAT=allocstat )
2516
         IF ( allocstat /= 0 ) THEN
2517
           CALL Fatal('GCR_Z','Failed to allocate memory of size: '&
2518
               //I2S(n)//' x '//I2S(m-1))
2519
         END IF
2520
         czero = CMPLX( 0.0_dp, 0.0_dp, KIND=dp )
2521
         V(1:n,1:m-1) = czero
2522
         S(1:n,1:m-1) = czero
2523
      END IF	
2524
      
2525
      CALL matvecsubr( x, r, ipar )
2526
      r(1:n) = b(1:n) - r(1:n)
2527
      
2528
      bnorm = normfun(n, b, 1)
2529
      rnorm = normfun(n, r, 1)
2530

2531
      IF (UseStopCFun) THEN
2532
        Residual = stopcfun(x,b,r,ipar,dpar)
2533
      ELSE
2534
        Residual = rnorm / bnorm
2535
      END IF
2536
      Converged = (Residual < MinTolerance) 
2537
      Diverged = (Residual > MaxTolerance) .OR. (Residual /= Residual)    
2538
      IF( Converged .OR. Diverged) RETURN
2539
      
2540
      DO k=1,Rounds
2541
	 !----------------------------------------------
2542
	 ! Check for restarting
2543
         !--------------------------------------------- 
2544
         IF ( MOD(k,m)==0 ) THEN
2545
            j = m
2546
         ELSE
2547
            j = MOD(k,m)
2548
            !--------------------------------------------
2549
            ! Compute the true residual when restarting:
2550
            !--------------------------------------------
2551
            IF ( (j==1) .AND. (k>1) ) THEN
2552
               CALL matvecsubr( x, r, ipar ) 
2553
               r(1:n) = b(1:n) - r(1:n)
2554
            END IF
2555
         END IF
2556
         !----------------------------------------------------------
2557
         ! Perform the preconditioning...
2558
         !---------------------------------------------------------------
2559
         CALL pcondrsubr( T1, r, ipar )         
2560
         CALL matvecsubr( T1, T2, ipar )
2561
         !--------------------------------------------------------------
2562
         ! Perform the orthogonalization of the search directions....
2563
         !--------------------------------------------------------------
2564
         DO i=1,j-1
2565
            beta = dotprodfun(n, V(1:n,i), 1, T2(1:n), 1 )
2566

2567
            T1(1:n) = T1(1:n) - beta * S(1:n,i)
2568
            T2(1:n) = T2(1:n) - beta * V(1:n,i)        
2569
         END DO
2570

2571
         alpha = normfun(n, T2(1:n), 1 )
2572
         T1(1:n) = 1.0d0/alpha * T1(1:n)
2573
         T2(1:n) = 1.0d0/alpha * T2(1:n)
2574

2575
         !-------------------------------------------------------------
2576
         ! The update of the solution and save the search data...
2577
         !------------------------------------------------------------- 
2578
         beta = dotprodfun(n, T2(1:n), 1, r(1:n), 1 )
2579
         x(1:n) = x(1:n) + beta * T1(1:n)      
2580
         r(1:n) = r(1:n) - beta * T2(1:n)
2581
	 IF ( j /= m ) THEN
2582
            S(1:n,j) = T1(1:n)
2583
            V(1:n,j) = T2(1:n)
2584
	 END IF       
2585

2586
         !--------------------------------------------------------------
2587
         ! Check whether the convergence criterion is met 
2588
         !--------------------------------------------------------------
2589
         rnorm = normfun(n, r, 1)
2590

2591
         IF (UseStopCFun) THEN
2592
           Residual = stopcfun(x,b,r,ipar,dpar)
2593
           IF( MOD(k,OutputInterval) == 0) THEN
2594
             WRITE (*, '(A, I6, 2E12.4)') '   gcr:',k, rnorm / bnorm, residual
2595
             CALL FLUSH(6)
2596
           END IF
2597
         ELSE
2598
           Residual = rnorm / bnorm
2599
           IF( MOD(k,OutputInterval) == 0) THEN
2600
             WRITE (*, '(A, I8, 3ES12.4,A)') '   gcrz:',k, residual, beta,'i'
2601
             CALL FLUSH(6)
2602
           END IF
2603
         END IF
2604
        
2605
         Converged = (Residual < MinTolerance)
2606
         !-----------------------------------------------------------------
2607
         ! Make an additional check that the true residual agrees with 
2608
         ! the iterated residual:
2609
         !-----------------------------------------------------------------
2610
         IF (Converged ) THEN
2611
            CALL matvecsubr( x, trueres, ipar )
2612
            trueres(1:n) = b(1:n) - trueres(1:n)
2613
            TrueResNorm = normfun(n, trueres, 1)
2614
            NormErr = ABS(TrueResNorm - rnorm)/TrueResNorm
2615
            IF ( NormErr > 1.0d-1 ) THEN
2616
               CALL Info('WARNING', 'Iterated GCR solution may not be accurate', Level=2)
2617
               WRITE( Message, * ) 'Iterated GCR residual norm = ', rnorm
2618
               CALL Info('WARNING', Message, Level=2)
2619
               WRITE( Message, * ) 'True residual norm = ', TrueResNorm
2620
               CALL Info('WARNING', Message, Level=2)   
2621
               CALL FLUSH(6)
2622
             END IF
2623
         END IF 
2624
         Diverged = (Residual > MaxTolerance) .OR. (Residual /= Residual)    
2625
         IF( Converged .OR. Diverged) EXIT
2626
        
2627
      END DO
2628
      
2629
      DEALLOCATE( R, T1, T2 )
2630
      IF ( m > 1 ) DEALLOCATE( S, V)
2631
      
2632
    END SUBROUTINE GCR_Z
2633
    
2634
!------------------------------------------------------------------------------
2635
  END SUBROUTINE itermethod_z_gcr
2636
!------------------------------------------------------------------------------
2637

2638

2639

2640
!------------------------------------------------------------------------------
2641
!> This routine provides a complex version of the BiCGstab(l) solver for linear systems.
2642
!------------------------------------------------------------------------------
2643
  SUBROUTINE itermethod_z_bicgstabl( xvec, rhsvec, &
2644
       ipar, dpar, work, matvecsubr, pcondlsubr, &
2645
       pcondrsubr, dotprodfun, normfun, stopcfun )
2646
 !------------------------------------------------------------------------------
2647
    USE huti_interfaces
2648
    IMPLICIT NONE
2649
    PROCEDURE( mv_iface_z ), POINTER :: matvecsubr
2650
    PROCEDURE( pc_iface_z ), POINTER :: pcondlsubr
2651
    PROCEDURE( pc_iface_z ), POINTER :: pcondrsubr
2652
    PROCEDURE( dotp_iface_z ), POINTER :: dotprodfun
2653
    PROCEDURE( norm_iface_z ), POINTER :: normfun
2654
    PROCEDURE( stopc_iface_z ), POINTER :: stopcfun
2655

2656
    INTEGER, DIMENSION(HUTI_IPAR_DFLTSIZE) :: ipar
2657
    COMPLEX(KIND=dp), TARGET, DIMENSION(HUTI_NDIM) :: xvec, rhsvec
2658
    DOUBLE PRECISION, DIMENSION(HUTI_DPAR_DFLTSIZE) :: dpar
2659
    COMPLEX(KIND=dp), DIMENSION(HUTI_WRKDIM,HUTI_NDIM) :: work
2660

2661
    COMPLEX(KIND=dp) :: y(HUTI_NDIM),f(HUTI_NDIM)
2662
    INTEGER :: ndim, i, PolynomialDegree
2663
    INTEGER :: Rounds, OutputInterval
2664
    REAL(KIND=dp) :: MinTol, MaxTol
2665
    LOGICAL :: Converged, Diverged
2666
    !--------------------------------------------------------------------------------    
2667

2668
    ndim = HUTI_NDIM
2669
    Rounds = HUTI_MAXIT
2670
    MinTol = HUTI_TOLERANCE
2671
    MaxTol = HUTI_MAXTOLERANCE
2672
    OutputInterval = HUTI_DBUGLVL
2673
    PolynomialDegree = HUTI_BICGSTABL_L 
2674
    !----------------------------------------------------------------------------
2675
    ! Transform the solution vector and the right-hand side vector to 
2676
    ! complex-valued vectors y and f
2677
    !---------------------------------------------------------------------------
2678
    DO i=1,ndim
2679
      y(i)=xvec(i)
2680
      f(i)=rhsvec(i)
2681
    END DO
2682

2683
    CALL ComplexBiCGStabl(ndim, GlobalMatrix, y, f, Rounds, MinTol, MaxTol, &
2684
         Converged, Diverged, OutputInterval, PolynomialDegree)
2685

2686
    IF(Converged) HUTI_INFO = HUTI_CONVERGENCE
2687
    IF(Diverged) HUTI_INFO = HUTI_DIVERGENCE
2688
    IF ( (.NOT. Converged) .AND. (.NOT. Diverged) ) HUTI_INFO = HUTI_MAXITER
2689

2690
    DO i=1,ndim
2691
      xvec(i) = y(i)
2692
    END DO
2693

2694
  CONTAINS 
2695

2696
    !-----------------------------------------------------------------------------------
2697
    SUBROUTINE ComplexBiCGStabl( n, A, x, b, MaxRounds, Tol, MaxTol, Converged, &
2698
         Diverged, OutputInterval, l)
2699
    !-----------------------------------------------------------------------------------
2700
    !   This subroutine solves complex linear systems Ax = b by using the BiCGStab(l) algorithm 
2701
    !   with l >= 2 and the right-oriented preconditioning. 
2702
    !
2703
    !   The subroutine has been written using as a starting point the work of D.R. Fokkema 
2704
    !   (subroutine zbistbl v1.1 1998). Dr. Fokkema has given the right to distribute
2705
    !   the derived work under GPL and hence the original copyright notice of the subroutine
2706
    !   has been removed accordingly.  
2707
    !
2708
    !----------------------------------------------------------------------------------- 
2709
      INTEGER :: l   ! polynomial degree
2710
      INTEGER :: n, MaxRounds, OutputInterval   
2711
      LOGICAL :: Converged, Diverged
2712
      TYPE(Matrix_t), POINTER :: A
2713
      COMPLEX(KIND=dp) :: x(n), b(n)
2714
      REAL(KIND=dp) :: Tol, MaxTol
2715
      !------------------------------------------------------------------------------
2716
      COMPLEX(KIND=dp) :: zzero, zone, t(n), kappa0, kappal 
2717
      REAL(KIND=dp) :: rnrm0, rnrm, mxnrmx, mxnrmr, errorind, &
2718
           delta = 1.0d-2, bnrm
2719
      INTEGER :: i, j, rr, r, u, xp, bp, z, zz, y0, yl, y, k, iwork(l-1), stat, Round
2720
      COMPLEX(KIND=dp) :: alpha, beta, omega, rho0, rho1, sigma, zdotc, varrho, hatgamma
2721
      COMPLEX(KIND=dp), ALLOCATABLE :: work(:,:), rwork(:,:)
2722
      LOGICAL rcmp, xpdt, EarlyExit
2723
      !------------------------------------------------------------------------------
2724

2725
      IF ( l < 2) CALL Fatal( 'RealBiCGStabl', 'Polynomial degree < 2' )
2726

2727
      IF ( ALL(x == CMPLX(0.0d0,0.0d0,kind=dp)) ) x = b
2728

2729
      zzero = CMPLX( 0.0d0,0.0d0, KIND=dp)
2730
      zone =  CMPLX( 1.0d0,0.0d0, KIND=dp)
2731

2732
      ALLOCATE( work(n,3+2*(l+1)), rwork(l+1,3+2*(l+1)) )
2733
      work = CMPLX( 0.0d0, 0.0d0, KIND=dp )
2734
      rwork = CMPLX( 0.0d0, 0.0d0, KIND=dp )
2735

2736
      rr = 1
2737
      r = rr+1
2738
      u = r+(l+1)
2739
      xp = u+(l+1)
2740
      bp = xp+1
2741

2742
      z = 1
2743
      zz = z+(l+1)
2744
      y0 = zz+(l+1)
2745
      yl = y0+1
2746
      y = yl+1
2747

2748
      CALL matvecsubr( x, work(1:n,r), ipar )
2749
      work(1:n,r) = b(1:n) - work(1:n,r)
2750
      bnrm = normfun(n, b(1:n), 1)
2751
      rnrm0 = normfun(n, work(1:n,r), 1)
2752

2753
      !-------------------------------------------------------------------
2754
      ! Check whether the initial guess satisfies the stopping criterion
2755
      !--------------------------------------------------------------------
2756
      errorind = rnrm0 / bnrm
2757
      Converged = (errorind < Tol)
2758
      Diverged = (errorind > MaxTol) .OR. (errorind /= errorind)
2759

2760
      IF( Converged .OR. Diverged) RETURN
2761
      EarlyExit = .FALSE.
2762

2763
      work(1:n,rr) = work(1:n,r) 
2764
      work(1:n,bp) = work(1:n,r)
2765
      work(1:n,xp) = x(1:n)
2766

2767
      rnrm = rnrm0
2768
      mxnrmx = rnrm0
2769
      mxnrmr = rnrm0
2770
      x(1:n) = zzero    
2771
      alpha = zzero
2772
      omega = zone
2773
      sigma = zone
2774
      rho0 = zone
2775

2776
      DO Round=1,MaxRounds 
2777
         !-------------------------
2778
         ! --- The BiCG part ---
2779
         !-------------------------
2780
         rho0 = -omega*rho0
2781

2782
         DO k=1,l
2783
            rho1 = dotprodfun(n, work(1:n,rr), 1, work(1:n,r+k-1), 1)
2784
            IF (rho0 == zzero) THEN
2785
               CALL Fatal( 'ComplexBiCGStab(l)', 'Breakdown error 1.' )
2786
            ENDIF
2787
            beta = alpha*(rho1/rho0)
2788
            rho0 = rho1
2789
            DO j=0,k-1
2790
               work(1:n,u+j) = work(1:n,r+j) - beta*work(1:n,u+j)
2791
            ENDDO
2792
            CALL pcondrsubr( t, work(1:n,u+k-1), ipar )
2793
            CALL matvecsubr( t, work(1:n,u+k),   ipar )
2794

2795
            sigma = dotprodfun(n, work(1:n,rr), 1, work(1:n,u+k), 1)
2796
            IF (sigma == zzero) THEN
2797
               CALL Fatal( 'ComplexBiCGStab(l)', 'Breakdown error 2.' )
2798
            ENDIF
2799
            alpha = rho1/sigma
2800
            x(1:n) = x(1:n) + alpha * work(1:n,u)
2801
            DO j=0,k-1
2802
               work(1:n,r+j) = work(1:n,r+j) - alpha * work(1:n,u+j+1)
2803
            ENDDO
2804
            CALL pcondrsubr( t, work(1:n,r+k-1), ipar )
2805
            CALL matvecsubr( t, work(1:n,r+k),   ipar )
2806
            rnrm = normfun(n, work(1:n,r), 1)
2807
            mxnrmx = MAX (mxnrmx, rnrm)
2808
            mxnrmr = MAX (mxnrmr, rnrm)
2809

2810
            !----------------------------------------------------------------------
2811
            ! In some simple cases, a few BiCG updates may already be enough to
2812
            ! obtain the solution. The following is for handling this special case. 
2813
            !----------------------------------------------------------------------
2814
            errorind = rnrm / bnrm
2815
            Converged = (errorind < Tol) 
2816
            IF (Converged) THEN
2817
               EarlyExit = .TRUE.
2818
               EXIT
2819
            END IF
2820

2821
         ENDDO
2822
         
2823
         IF (EarlyExit) EXIT        
2824

2825
         !--------------------------------------
2826
         ! --- The convex polynomial part ---
2827
         !--------------------------------------
2828

2829
         DO i=1,l+1
2830
            DO j=1,i
2831
               rwork(i,j) = dotprodfun(n, work(1:n,r+i-1), 1, work(1:n,r+j-1),1 ) 
2832
            END DO
2833
         END DO
2834
         DO j=2,l+1
2835
            rwork(1:j-1,j) = CONJG( rwork(j,1:j-1) )
2836
         END DO
2837

2838
         rwork(1:l+1,zz:zz+l) = rwork(1:l+1,z:z+l)
2839
         CALL zgetrf (l-1, l-1, rwork(2:l,zz+1:zz+l-1), l-1, &
2840
              iwork, stat)
2841

2842
         ! --- tilde r0 and tilde rl (small vectors)
2843

2844
         rwork(1,y0) = -zone
2845
         rwork(2:l,y0) = rwork(2:l,z) 
2846
         CALL zgetrs('n', l-1, 1, rwork(2:l,zz+1:zz+l-1), l-1, iwork, &
2847
              rwork(2:l,y0), l-1, stat)
2848
         rwork(l+1,y0) = zzero
2849

2850
         rwork(1,yl) = zzero
2851
         rwork(2:l,yl) = rwork(2:l,z+l) 
2852
         CALL zgetrs ('n', l-1, 1, rwork(2:l,zz+1:zz+l-1), l-1, iwork, &
2853
              rwork(2:l,yl), l-1, stat)
2854
         rwork(l+1,yl) = -zone
2855

2856
         ! --- Convex combination
2857
         call zmv( rwork(1:l+1,z:z+l), rwork(1:l+1,y0), rwork(1:l+1,y), l+1 )
2858
         kappa0 = SQRT( ABS(zdotc(l+1, rwork(1:l+1,y0), 1, rwork(1:l+1,y), 1)) )  ! replace zdotc
2859

2860
         call zmv( rwork(1:l+1,z:z+l), rwork(1:l+1,yl), rwork(1:l+1,y), l+1 )
2861
         kappal = SQRT( ABS(zdotc(l+1, rwork(1:l+1,yl), 1, rwork(1:l+1,y), 1)) )  ! replace zdotc
2862

2863
         call zmv( rwork(1:l+1,z:z+l), rwork(1:l+1,y0), rwork(1:l+1,y), l+1 )
2864
         varrho = zdotc(l+1, rwork(1:l+1,yl), 1, rwork(1:l+1,y), 1) / &           ! replace zdotc
2865
              (kappa0*kappal)
2866

2867
         hatgamma = varrho/ABS(varrho) * MAX(ABS(varrho),7d-1) * &
2868
              kappa0/kappal
2869
         rwork(1:l+1,y0) = rwork(1:l+1,y0) - hatgamma * rwork(1:l+1,yl)
2870

2871
         !  --- Update
2872

2873
         omega = rwork(l+1,y0)
2874
         DO j=1,l
2875
            work(1:n,u) = work(1:n,u) - rwork(j+1,y0) * work(1:n,u+j)
2876
            x(1:n) = x(1:n) + rwork(j+1,y0) * work(1:n,r+j-1)
2877
            work(1:n,r) = work(1:n,r) - rwork(j+1,y0) * work(1:n,r+j)
2878
         ENDDO
2879

2880
         call zmv( rwork(1:l+1,z:z+l), rwork(1:l+1,y0), rwork(1:l+1,y), l+1 )
2881
         rnrm = SQRT( ABS(zdotc(l+1, rwork(1:l+1,y0), 1, rwork(1:l+1,y), 1)) )
2882

2883
         !---------------------------------------
2884
         !  --- The reliable update part ---
2885
         !---------------------------------------
2886

2887
         mxnrmx = MAX (mxnrmx, rnrm)
2888
         mxnrmr = MAX (mxnrmr, rnrm)
2889
         xpdt = (rnrm < delta*rnrm0 .AND. rnrm0 < mxnrmx)
2890
         rcmp = ((rnrm < delta*mxnrmr .AND. rnrm0 < mxnrmr) .OR. xpdt)
2891
         IF (rcmp) THEN
2892
            ! PRINT *, 'Performing residual update...'
2893
            CALL pcondrsubr( t, x, ipar )
2894
            CALL matvecsubr( t, work(1:n,r), ipar )
2895
            work(1:n,r) = work(1:n,bp) - work(1:n,r)
2896
            mxnrmr = rnrm
2897
            IF (xpdt) THEN
2898
               ! PRINT *, 'Performing solution update...'
2899
               work(1:n,xp) = work(1:n,xp) + t(1:n)
2900
               x(1:n) = zzero
2901
               work(1:n,bp) = work(1:n,r)
2902
               mxnrmx = rnrm
2903
            ENDIF
2904
         ENDIF
2905

2906
         IF (rcmp) THEN
2907
            IF (xpdt) THEN       
2908
               t(1:n) = work(1:n,xp)
2909
            ELSE
2910
               t(1:n) = t(1:n) + work(1:n,xp)  
2911
            END IF
2912
         ELSE
2913
            CALL pcondrsubr( t, x, ipar )
2914
            t(1:n) = t(1:n)+work(1:n,xp)
2915
         END IF
2916

2917
         errorind = rnrm/bnrm
2918
         IF( MOD(Round,OutputInterval) == 0) THEN
2919
           WRITE (*, '(I8, E11.4)') Round, errorind
2920
           CALL FLUSH(6)
2921
         END IF
2922

2923
         Converged = (errorind < Tol) 
2924
         Diverged = (errorind > MaxTol) .OR. (errorind /= errorind)    
2925
         IF( Converged .OR. Diverged) EXIT    
2926
      END DO
2927

2928
      IF( EarlyExit .AND. (OutputInterval/=HUGE(OutputInterval)) ) THEN
2929
        WRITE (*, '(I8, E11.4)') Round, errorind
2930
        CALL FLUSH(6)
2931
      END IF
2932

2933
      !------------------------------------------------------------
2934
      ! We have solved z = P*x, so finally solve the true unknown x
2935
      !------------------------------------------------------------
2936
      t(1:n) = x(1:n)
2937
      CALL pcondrsubr( x, t, ipar )
2938
      x(1:n) = x(1:n) + work(1:n,xp)      
2939

2940
    !----------------------------------------------------------
2941
    END SUBROUTINE ComplexBiCGStabl
2942
    !----------------------------------------------------------
2943

2944
!--------------------------------------------------------------
2945

2946
    SUBROUTINE zmv(a,x,y,n)
2947
      integer, INTENT(in) :: n
2948
      complex(kind=dp), INTENT(in)  ::a(n,n), x(n)
2949
      complex(kind=dp), INTENT(out) ::y(n)
2950
 
2951
      complex(kind=dp), parameter :: zone  = cmplx(1._dp, 0._dp,kind=dp)
2952
      complex(kind=dp), parameter :: zzero = cmplx(0._dp, 0._dp,kind=dp)
2953
    
2954
      IF(n>8) THEN
2955
        CALL zhemv ('u', n, zone, a, n, x, 1, zzero, y, 1)
2956
        RETURN
2957
      END IF
2958

2959
      y(1) = a(1,1)*x(1) + a(1,2)*x(2) + a(1,3)*x(3)
2960
      y(2) = a(2,1)*x(1) + a(2,2)*x(2) + a(2,3)*x(3)
2961
      y(3) = a(3,1)*x(1) + a(3,2)*x(2) + a(3,3)*x(3)
2962
      IF(n==3) RETURN
2963

2964
      y(1) = y(1) + a(1,4)*x(4)
2965
      y(2) = y(2) + a(2,4)*x(4)
2966
      y(3) = y(3) + a(3,4)*x(4)
2967
      y(4) = a(4,1)*x(1)+a(4,2)*x(2)+a(4,3)*x(3)+a(4,4)*x(4)
2968
      IF(n==4) RETURN
2969

2970
      y(1) = y(1) + a(1,5)*x(5)
2971
      y(2) = y(2) + a(2,5)*x(5)
2972
      y(3) = y(3) + a(3,5)*x(5)
2973
      y(4) = y(4) + a(4,5)*x(5)
2974
      y(5) = a(5,1)*x(1)+a(5,2)*x(2)+a(5,3)*x(3)+a(5,4)*x(4)+ &
2975
             a(5,5)*x(5)
2976
      IF(n==5) RETURN
2977

2978
      y(1) = y(1) + a(1,6)*x(6)
2979
      y(2) = y(2) + a(2,6)*x(6)
2980
      y(3) = y(3) + a(3,6)*x(6)
2981
      y(4) = y(4) + a(4,6)*x(6)
2982
      y(5) = y(5) + a(5,6)*x(6)
2983
      y(6) = a(6,1)*x(1)+a(6,2)*x(2)+a(6,3)*x(3)+a(6,4)*x(4)+ &
2984
             a(6,5)*x(5)+a(6,6)*x(6)
2985
      IF(n==6) RETURN
2986

2987
      y(1) = y(1) + a(1,7)*x(7)
2988
      y(2) = y(2) + a(2,7)*x(7)
2989
      y(3) = y(3) + a(3,7)*x(7)
2990
      y(4) = y(4) + a(4,7)*x(7)
2991
      y(5) = y(5) + a(5,7)*x(7)
2992
      y(6) = y(6) + a(6,7)*x(7)
2993
      y(7) = a(7,1)*x(1)+a(7,2)*x(2)+a(7,3)*x(3)+a(7,4)*x(4)+ &
2994
             a(7,5)*x(5)+a(7,6)*x(6)+a(7,7)*x(7)
2995
      IF(n==7) RETURN
2996

2997
      y(1) = y(1) + a(1,8)*x(8)
2998
      y(2) = y(2) + a(2,8)*x(8)
2999
      y(3) = y(3) + a(3,8)*x(8)
3000
      y(4) = y(4) + a(4,8)*x(8)
3001
      y(5) = y(5) + a(5,8)*x(8)
3002
      y(6) = y(6) + a(6,8)*x(8)
3003
      y(7) = y(7) + a(7,8)*x(8)
3004
      y(8) = a(8,1)*x(1)+a(8,2)*x(2)+a(8,3)*x(3)+a(8,4)*x(4)+ &
3005
             a(8,5)*x(5)+a(8,6)*x(6)+a(8,7)*x(7)+a(8,8)*x(8)
3006

3007
    END SUBROUTINE zmv
3008

3009
!--------------------------------------------------------------
3010
  END SUBROUTINE itermethod_z_bicgstabl
3011
!--------------------------------------------------------------
3012

3013

3014
!------------------------------------------------------------------------------
3015
!> This routine provides a complex version of the IDR(s) solver for linear systems.
3016
!------------------------------------------------------------------------------
3017
  SUBROUTINE itermethod_z_idrs( xvec, rhsvec, &
3018
      ipar, dpar, work, matvecsubr, pcondlsubr, &
3019
      pcondrsubr, dotprodfun, normfun, stopcfun )
3020
!------------------------------------------------------------------------------
3021
    USE huti_interfaces
3022
    IMPLICIT NONE
3023
    PROCEDURE( mv_iface_z ), POINTER :: matvecsubr
3024
    PROCEDURE( pc_iface_z ), POINTER :: pcondlsubr
3025
    PROCEDURE( pc_iface_z ), POINTER :: pcondrsubr
3026
    PROCEDURE( dotp_iface_z ), POINTER :: dotprodfun
3027
    PROCEDURE( norm_iface_z ), POINTER :: normfun
3028
    PROCEDURE( stopc_iface_z ), POINTER :: stopcfun
3029

3030
    INTEGER, DIMENSION(HUTI_IPAR_DFLTSIZE) :: ipar
3031
    COMPLEX(KIND=dp), TARGET, DIMENSION(HUTI_NDIM) :: xvec, rhsvec
3032
    DOUBLE PRECISION, DIMENSION(HUTI_DPAR_DFLTSIZE) :: dpar
3033
    COMPLEX(KIND=dp), DIMENSION(HUTI_WRKDIM,HUTI_NDIM) :: work
3034

3035
    COMPLEX(KIND=dp) :: y(HUTI_NDIM),f(HUTI_NDIM)
3036
    INTEGER :: ndim, i, s
3037
    INTEGER :: Rounds, OutputInterval
3038
    REAL(KIND=dp) :: MinTol, MaxTol
3039
    LOGICAL :: Converged, Diverged
3040
    !--------------------------------------------------------------------------------    
3041

3042
    ndim = HUTI_NDIM
3043
    Rounds = HUTI_MAXIT
3044
    MinTol = HUTI_TOLERANCE
3045
    MaxTol = HUTI_MAXTOLERANCE
3046
    OutputInterval = HUTI_DBUGLVL
3047
    s = HUTI_IDRS_S 
3048
    !----------------------------------------------------------------------------
3049
    ! Transform the solution vector and the right-hand side vector to 
3050
    ! complex-valued vectors y and f
3051
    !---------------------------------------------------------------------------
3052
    DO i=1,ndim
3053
      y(i)=xvec(i)
3054
      f(i)=rhsvec(i)
3055
    END DO
3056
    CALL ComplexIDRS(ndim, GlobalMatrix, y, f, Rounds, MinTol, MaxTol, &
3057
         Converged, Diverged, OutputInterval, s )
3058

3059
    IF(Converged) HUTI_INFO = HUTI_CONVERGENCE
3060
    IF(Diverged) HUTI_INFO = HUTI_DIVERGENCE
3061
    IF ( (.NOT. Converged) .AND. (.NOT. Diverged) ) HUTI_INFO = HUTI_MAXITER
3062

3063
    DO i=1,ndim
3064
      xvec(i) = y(i)
3065
    END DO
3066

3067
  CONTAINS 
3068

3069
!----------------------------------------------------------------------------------- 
3070
!   This subroutine solves complex linear systems Ax = b by using the IDR(s) algorithm 
3071
!   with s >= 1 and the right-oriented preconditioning. 
3072
!
3073
!   The subroutine ComplexIDRS has been written by M.B. van Gijzen
3074
!-----------------------------------------------------------------------------------
3075
    SUBROUTINE ComplexIDRS( n, A, x, b, MaxRounds, Tol, MaxTol, Converged, &
3076
        Diverged, OutputInterval, s )
3077
!-----------------------------------------------------------------------------------
3078
      INTEGER :: s  
3079
      INTEGER :: n, MaxRounds, OutputInterval   
3080
      LOGICAL :: Converged, Diverged, UseStopCFun
3081
      TYPE(Matrix_t), POINTER :: A
3082
      COMPLEX(KIND=dp) :: x(n), b(n)
3083
      REAL(KIND=dp) :: Tol, MaxTol
3084
!------------------------------------------------------------------------------
3085

3086
      ! Local arrays:
3087
!     REAL(kind=dp) :: Pr(n,s), Pi(n,s) 
3088
!     COMPLEX(kind=dp) :: P(n,s)
3089
!     COMPLEX(kind=dp) :: G(n,s)
3090
!     COMPLEX(kind=dp) :: U(n,s)
3091
!     COMPLEX(kind=dp) :: r(n) 
3092
!     COMPLEX(kind=dp) :: v(n)   
3093
!     COMPLEX(kind=dp) :: t(n)  
3094
!     COMPLEX(kind=dp) :: M(s,s), f(s), mu(s)
3095
!     COMPLEX(kind=dp) :: alpha(s), beta(s), gamma(s)
3096

3097
      REAL(kind=dp), ALLOCATABLE :: Pr(:,:), Pi(:,:) 
3098
      COMPLEX(kind=dp), ALLOCATABLE :: P(:,:)
3099
      COMPLEX(kind=dp), ALLOCATABLE :: G(:,:)
3100
      COMPLEX(kind=dp), ALLOCATABLE :: U(:,:)
3101
      COMPLEX(kind=dp), ALLOCATABLE :: r(:) 
3102
      COMPLEX(kind=dp), ALLOCATABLE :: v(:)   
3103
      COMPLEX(kind=dp), ALLOCATABLE :: t(:)  
3104
      COMPLEX(kind=dp), ALLOCATABLE :: M(:,:), f(:), mu(:)
3105
      COMPLEX(kind=dp), ALLOCATABLE :: alpha(:), beta(:), gamma(:)
3106

3107
      COMPLEX(kind=dp) :: om, tr    
3108
      REAL(kind=dp) :: nr, nt, rho, kappa
3109

3110
      INTEGER :: iter                         ! number of iterations
3111
      INTEGER :: ii                           ! inner iterations index
3112
      INTEGER :: jj                           ! G-space index
3113
      REAL(kind=dp) :: normb, normr, errorind ! for tolerance check
3114
      INTEGER :: i,j,k,l                      ! loop counters
3115

3116
      UseStopCFun = HUTI_STOPC == HUTI_USUPPLIED_STOPC
3117
      
3118
      ALLOCATE( Pr(n,s), Pi(n,s), P(n,s), G(n,s), U(n,s), r(n), v(n), t(n), &
3119
            M(s,s), f(s), mu(s), alpha(s), beta(s), gamma(s))
3120
      U = 0.0d0
3121

3122
      ! Compute initial residual, set absolute tolerance
3123
      normb = normfun(n,b,1)
3124
      CALL matvecsubr( x, t, ipar )
3125
      r = b - t
3126
      IF (UseStopCFun) THEN
3127
        errorind = stopcfun(x,b,r,ipar,dpar)
3128
      ELSE
3129
        normr = normfun(n,r,1)
3130
        errorind = normr / normb
3131
      END IF
3132
      
3133
      !-------------------------------------------------------------------
3134
      ! Check whether the initial guess satisfies the stopping criterion
3135
      !--------------------------------------------------------------------
3136
      Converged = (errorind < Tol)
3137
      Diverged = (errorind > MaxTol) .OR. (errorind /= errorind)
3138

3139
      IF ( Converged .OR. Diverged ) RETURN
3140

3141
      ! Define P and kappa 
3142
      CALL RANDOM_NUMBER(Pr)
3143
      CALL RANDOM_NUMBER(Pi)
3144
      P = Pr + (0.,1.)*Pi
3145

3146
      DO j = 1,s
3147
        DO k = 1,j-1
3148
          alpha(k) = dotprodfun(n, P(:,k), 1, P(:,j), 1 )
3149
          P(:,j) = P(:,j) - alpha(k)*P(:,k)
3150
        END DO
3151
        P(:,j) = P(:,j)/normfun(n,P(:,j),1)
3152
      END DO
3153
      kappa = 0.7d0
3154

3155
      ! Initialize local variables:
3156
      M = 0.0d0
3157
      om = 1.0d0
3158
      iter = 0
3159
      jj = 0
3160
      ii = 0
3161
      ! This concludes the initialisation phase
3162

3163
      ! Main iteration loop, build G-spaces:
3164

3165
      DO WHILE ( .NOT. Converged .AND. .NOT. Diverged )  ! start of iteration loop
3166

3167
        !!+++++++++++++++++++++++++++++++++++++++++++++++++++++++
3168
        ! Generate s vectors in G_j
3169
        !!+++++++++++++++++++++++++++++++++++++++++++++++++++++++
3170

3171
        ! New right-hand side for small system:
3172
        DO k = 1,s
3173
          f(k) = dotprodfun(n, P(:,k), 1, r, 1 )
3174
        END DO
3175

3176
        DO k = 1,s
3177

3178
          ! Update inner iteration counter
3179
          ii = ii + 1
3180

3181
          ! Compute new v
3182
          v = r 
3183
          IF ( jj > 0 ) THEN
3184

3185
            ! Solve small system (Note: M is lower triangular) and make v orthogonal to P:
3186
            DO i = k,s
3187
              gamma(i) = f(i)
3188
              DO j = k,i-1
3189
                gamma(i) = gamma(i) - M(i,j)*gamma(j)
3190
              END DO
3191
              gamma(i) = gamma(i)/M(i,i)
3192
              v = v - gamma(i)*G(:,i)
3193
            END DO
3194

3195
            ! Compute new U(:,k)
3196
            CALL pcondrsubr( t, v, ipar )
3197
            t = om*t
3198
            DO i = k,s
3199
              t = t + gamma(i)*U(:,i)
3200
            END DO
3201
            U(:,k) = t
3202

3203
          ELSE 
3204

3205
            ! Updates for the first s iterations (in G_0):
3206
            CALL pcondrsubr( U(:,k), v, ipar )
3207

3208
          END IF
3209

3210
          ! Compute new G(:,k), G(:,k) is in space G_j
3211
          CALL matvecsubr( U(:,k), G(:,k), ipar )
3212

3213
          ! Bi-Orthogonalise the new basis vectors: 
3214
          DO i = 1,s
3215
            mu(i) = dotprodfun(n, P(:,i), 1, G(:,k), 1 )
3216
          END DO
3217
          DO i = 1,k-1
3218
            alpha(i) = mu(i)
3219
            DO j = 1, i-1
3220
              alpha(i) = alpha(i) - M(i,j)*alpha(j)
3221
            END DO
3222
            alpha(i) = alpha(i)/M(i,i)
3223
            G(:,k) = G(:,k) - G(:,i)*alpha(i)
3224
            U(:,k) = U(:,k) - U(:,i)*alpha(i)
3225
            mu(k:s)  = mu(k:s)  - M(k:s,i)*alpha(i)
3226
          END DO
3227
          M(k:s,k) = mu(k:s)
3228

3229
          ! Break down?
3230
          IF ( ABS(M(k,k)) <= TINY(tol) ) THEN
3231
            Diverged = .TRUE.
3232
            EXIT
3233
          END IF
3234

3235
          ! Make r orthogonal to p_i, i = 1..k, update solution and residual 
3236
          beta(k) = f(k)/M(k,k)
3237
          r = r - beta(k)*G(:,k)
3238
          x = x + beta(k)*U(:,k)
3239

3240
          ! New f = P'*r (first k  components are zero)
3241
          IF ( k < s ) THEN
3242
            f(k+1:s)   = f(k+1:s) - beta(k)*M(k+1:s,k)
3243
          END IF
3244

3245
          ! Check for convergence
3246
          IF (UseStopCFun) THEN
3247
            errorind = stopcfun(x,b,r,ipar,dpar)
3248
          ELSE
3249
            normr = normfun(n,r,1)
3250
            errorind = normr/normb
3251
          END IF
3252
          iter = iter + 1
3253
          
3254
          IF( MOD(iter,OutputInterval) == 0) THEN
3255
            WRITE (*, '(I8, E11.4)') iter, errorind
3256
            CALL FLUSH(6)
3257
          END IF
3258

3259
          Converged = (errorind < Tol)
3260
          Diverged = (errorind > MaxTol) .OR. (errorind /= errorind)
3261
          IF ( Converged .OR. Diverged ) EXIT
3262
          IF (iter == MaxRounds) EXIT
3263

3264
        END DO ! Now we have computed s+1 vectors in G_j
3265
        IF ( Converged .OR. Diverged ) EXIT
3266
        IF (iter == MaxRounds) EXIT
3267

3268
        !!+++++++++++++++++++++++++++++++++++++++++++++++++++++++
3269
        ! Compute first residual in G_j+1
3270
        !!+++++++++++++++++++++++++++++++++++++++++++++++++++++++
3271

3272
        ! Update G-space counter
3273
        jj = jj + 1
3274

3275
        ! Compute first residual in G_j+1
3276
        ! Note: r is already perpendicular to P so v = r
3277

3278
        ! Preconditioning:
3279
        CALL pcondrsubr( v, r, ipar )
3280
        ! Matrix-vector multiplication:
3281
        CALL matvecsubr( v, t, ipar )
3282

3283
        ! Computation of a new omega
3284
        ! 'Maintaining the convergence':
3285
        nr = normfun(n,r,1)
3286
        nt = normfun(n,t,1)
3287
        tr = dotprodfun(n, t, 1, r, 1 )
3288
        rho = ABS(tr/(nt*nr))
3289
        om=tr/(nt*nt)
3290
        IF ( rho < kappa ) THEN
3291
          om = om*kappa/rho
3292
        END IF
3293

3294
        IF ( ABS(om) <= EPSILON(tol) ) THEN 
3295
          Diverged = .TRUE.
3296
          EXIT
3297
        END IF
3298

3299
        ! Update solution and residual
3300
        r = r - om*t 
3301
        x = x + om*v 
3302

3303
        ! Check for convergence
3304
        IF (UseStopCFun) THEN
3305
          errorind = stopcfun(x,b,r,ipar,dpar)
3306
        ELSE
3307
          normr = normfun(n,r,1)
3308
          errorind = normr/normb
3309
        END IF
3310
        iter = iter + 1
3311

3312
        IF( MOD(iter,OutputInterval) == 0) THEN
3313
          WRITE (*, '(I8, E11.4)') iter, errorind
3314
          CALL FLUSH(6)
3315
        END IF
3316

3317
        Converged = (errorind < Tol)
3318
        Diverged = (errorind > MaxTol) .OR. (errorind /= errorind)
3319
        IF (iter == MaxRounds) EXIT
3320

3321
      END DO ! end of while loop
3322

3323
    !----------------------------------------------------------
3324
    END SUBROUTINE ComplexIDRS
3325
    !----------------------------------------------------------
3326

3327
!--------------------------------------------------------------
3328
  END SUBROUTINE itermethod_z_idrs
3329
!--------------------------------------------------------------
3330

3331
END MODULE IterativeMethods
3332

3333
!> \} ElmerLib
3334

3335
Product

Resources

Company
