CoCalc -- huti_gmres.F90

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

27
  !
28
  ! Subroutines to implement Generalized Minimum Residual iterative method
29
  !
30

31
#include  "huti_fdefs.h" 
32

33
  !*************************************************************************
34
  !*************************************************************************
35
  !
36
  ! These subroutines are based on a book by Barret et al.:
37
  ! "Templates for the Solution of Linear Systems: Building Blocks for
38
  !  Iterative Methods", 1993.
39
  !
40
  ! All matrix-vector operations are done externally, so we do not need
41
  ! to know about the matrix structure (sparse or dense). So has the
42
  ! memory allocation for the working arrays done also externally.
43

44
  !*************************************************************************
45
  ! Work array is used in the following order:
46
  ! work(:,1) = z
47
  ! work(:,2) = p
48
  ! work(:,3) = q
49
  ! work(:,4) = r
50
  !
51
  ! Definitions to make the code more understandable and to make it look
52
  ! like the pseudo code (these are common to all precisions)
53
  !
54

55
#define  X  xvec 
56
#define  B  rhsvec 
57

58
#define  W  work(:,1) 
59
#define  W_ind  1 
60
#define  R  work(:,2) 
61
#define  R_ind  2 
62
#define  S  work(:,3) 
63
#define  S_ind  3 
64
#define  VTMP  work(:,4) 
65
#define  VTMP_ind  4 
66
#define  T1V  work(:,5) 
67
#define  T1V_ind  5 
68
#define  V  work(:,6) 
69
#define  V_ind  6 
70

71
contains
72

73
  !*************************************************************************
74

75
  !*************************************************************************
76
  !*************************************************************************
77
  ! Single precision version
78
  !*************************************************************************
79
  !*************************************************************************
80

81
  subroutine  huti_sgmressolv  ( ndim, wrkdim, xvec, rhsvec, &
82
       ipar, dpar, work, matvecsubr, pcondlsubr, &
83
       pcondrsubr, dotprodfun, normfun, stopcfun )
84

85
    use huti_interfaces
86
    implicit none
87

88
    procedure( mv_iface_s ), pointer :: matvecsubr
89
    procedure( pc_iface_s ), pointer :: pcondlsubr
90
    procedure( pc_iface_s ), pointer :: pcondrsubr
91
    procedure( dotp_iface_s ), pointer :: dotprodfun
92
    procedure( norm_iface_s ), pointer :: normfun
93
    procedure( stopc_iface_s ), pointer :: stopcfun
94

95
    ! Parameters
96

97
    integer :: ndim, wrkdim
98
    real, dimension(ndim) :: xvec, rhsvec
99
    integer, dimension(HUTI_IPAR_DFLTSIZE) :: ipar
100
    double precision, dimension(HUTI_DPAR_DFLTSIZE) :: dpar
101
    real, dimension(ndim,wrkdim) :: work
102

103
    ! Local variables
104

105
    integer :: iter_count
106
    real :: residual, rhsnorm, precrhsnorm
107
    real :: bnrm, alpha, beta
108

109
    INTEGER :: reit, info, i, j, k, l, m, n
110

111
    real :: temp, temp2, error, gamma
112

113
    ! Local arrays
114

115
    real, DIMENSION(HUTI_GMRES_RESTART+1,HUTI_GMRES_RESTART+1) :: H
116
    real, &
117
         DIMENSION((HUTI_GMRES_RESTART+1)*(HUTI_GMRES_RESTART+1)) :: HLU
118
    real, DIMENSION(HUTI_GMRES_RESTART+1) :: CS, SN, Y
119

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

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

131
    iter_count = 1
132

133
    bnrm = normfun( HUTI_NDIM, B, 1 )
134

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

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

146
    ! The following applies for all matrix operations in this solver
147

148
    HUTI_EXTOP_MATTYPE = HUTI_MAT_NOTTRPSED
149

150
    ! Generate vector X if needed
151

152
    if ( HUTI_INITIALX .eq. HUTI_RANDOMX ) then
153
       call  huti_srandvec   ( X, ipar )
154
    else if ( HUTI_INITIALX .ne. HUTI_USERSUPPLIEDX ) then
155
       X = 1
156
    end if
157

158
    call pcondrsubr( T1V, X, ipar )
159
    call matvecsubr( T1V, R, ipar )
160
    T1V = B - R
161
    call pcondlsubr( R, T1V, ipar )
162

163
    m = HUTI_GMRES_RESTART
164
    work(:,V_ind+1-1:V_ind+m+1-1) = 0
165
    H = 0
166
    CS = 0
167
    SN = 0
168
    VTMP = 0
169
    work(1,VTMP_ind) = 1.0
170

171
    !
172
    ! This is where the loop starts (that is we continue from here after
173
    ! the first iteration)
174
    !
175

176
300 continue
177

178
    call pcondrsubr( T1V, X, ipar )
179
    call matvecsubr( T1V, R, ipar )
180
    T1V = B - R
181
    call pcondlsubr( R, T1V, ipar )
182

183
    alpha = normfun( HUTI_NDIM, R, 1 )
184
    if ( alpha .eq. 0 ) then
185
       HUTI_INFO = HUTI_GMRES_ALPHA
186
       go to 1000
187
    end if
188

189
    work(:,V_ind+1-1) = R / alpha
190
    S = alpha * VTMP
191

192
    !
193
    ! Construct orthonormal
194
    !
195

196
    DO i = 1, m
197
       call pcondrsubr( W, work(:,V_ind+i-1), ipar )
198
       call matvecsubr( W, T1V, ipar )
199
       call pcondlsubr( W, T1V, ipar )
200

201
       DO k = 1, i
202
          H(k,i) = dotprodfun( HUTI_NDIM, W, 1, work(:,V_ind+k-1), 1 )
203
          W = W - H(k,i) * work(:,V_ind+k-1)
204
       END DO
205

206
       beta = normfun( HUTI_NDIM, W, 1 )
207
       if ( beta .eq. 0 ) then
208
          HUTI_INFO = HUTI_GMRES_BETA
209
          go to 1000
210
       end if
211

212
       H(i+1,i) = beta
213
       work(:,V_ind+i+1-1) = W / H(i+1, i)
214

215
       !
216
       ! Compute the Givens rotation
217
       !
218

219
       DO k = 1, i-1
220
          temp = CS(k) * H(k,i) + SN(k) * H(k+1,i)
221
          H(k+1,i) = -1 * SN(k) * H(k,i) + CS(k) * H(k+1,i)
222
          H(k,i) = temp
223
       END DO
224

225
       IF ( H(i+1,i) .eq. 0 ) THEN
226
          CS(i) = 1; SN(i) = 0
227
       ELSE
228
          IF ( abs( H(i+1,i) ) .gt. abs( H(i,i) ) ) THEN
229
             temp2 = H(i,i) / H(i+1,i)
230
             SN(i) = 1 / sqrt( 1 + ( temp2 * temp2 ))
231
             CS(i) = temp2 * SN(i)
232
          ELSE
233
             temp2 = H(i+1,i) / H(i,i)
234
             CS(i) = 1 / sqrt( 1 + ( temp2 * temp2 ))
235
             SN(i) = temp2 * CS(i)
236
          END IF
237
       END IF
238

239
       temp = CS(i) * work(i,S_ind)
240
       work(i+1,S_ind) = -1 * SN(i) * work(i,S_ind)
241
       work(i,S_ind) = temp
242

243
       H(i,i) = ( CS(i) * H(i,i) ) + ( SN(i) * H(i+1,i) )
244
       H(i+1,i) = 0
245

246
       error = abs( work(i+1,S_ind) ) / bnrm 
247
       IF ( REAL( error ) .lt.  HUTI_TOLERANCE ) THEN
248

249
          HLU = 0; j = 1
250
          do k = 1, i
251
             do l = 1, i
252
                HLU(j) = H(l,k)
253
                j = j + 1
254
             end do
255
          end do
256

257
          call  huti_slusolve  ( i, HLU, Y, work(:,S_ind) )
258

259
          X = X + MATMUL( work(:,V_ind+1-1:V_ind+i-1), Y(1:i) )
260

261
          EXIT
262
       END IF
263

264
    END DO
265

266
    IF ( REAL( error ) .lt. HUTI_TOLERANCE ) THEN
267
       GOTO 500
268
    END IF
269

270
    HLU = 0; j = 1
271
    do k = 1, m
272
       do l = 1, m
273
          HLU(j) = H(l,k)
274
          j = j + 1
275
       end do
276
    end do
277

278
    call  huti_slusolve  ( m, HLU, Y, work(:,S_ind) )
279

280
    X = X + MATMUL( work(:,V_ind+1-1:V_ind+m-1), Y(1:m) )
281

282
500 CONTINUE
283

284
    !
285
    ! Check the convergence
286
    !
287

288
    select case (HUTI_STOPC)
289
    case (HUTI_TRUERESIDUAL)
290
       call pcondrsubr( T1V, X, ipar )
291
       call matvecsubr( T1V, R, ipar )
292
       T1V = B - R
293
       call pcondlsubr( R, T1V, ipar )
294
       residual = normfun( HUTI_NDIM, R, 1 )
295
    case (HUTI_TRESID_SCALED_BYB)
296
       call pcondrsubr( T1V, X, ipar )
297
       call matvecsubr( T1V, R, ipar )
298
       T1V = B - R
299
       call pcondlsubr( R, T1V, ipar )
300
       residual = normfun( HUTI_NDIM, R, 1 ) / rhsnorm
301
    case (HUTI_PSEUDORESIDUAL)
302
       call matvecsubr( X, R, ipar )
303
       R = R - B
304
       call pcondlsubr( T1V, R, ipar )
305
       residual = normfun( HUTI_NDIM, T1V, 1 )
306
    case (HUTI_PRESID_SCALED_BYB)
307
       call matvecsubr( X, R, ipar )
308
       R = R - B
309
       call pcondlsubr( T1V, R, ipar )
310
       residual = normfun( HUTI_NDIM, T1V, 1 ) / rhsnorm
311
    case (HUTI_PRESID_SCALED_BYPRECB)
312
       call matvecsubr( X, R, ipar )
313
       R = R - B
314
       call pcondlsubr( T1V, R, ipar )
315
       residual = normfun( HUTI_NDIM, T1V, 1 ) / precrhsnorm
316
    case (HUTI_XDIFF_NORM)
317
       T1V = MATMUL( work(:,V_ind+1-1:V_ind+m-1), Y(1:m) )
318
       residual = normfun( HUTI_NDIM, T1V, 1 )
319
    case (HUTI_USUPPLIED_STOPC)
320
       residual = stopcfun( X, B, R, ipar, dpar )
321
    case default
322
       call pcondrsubr( T1V, X, ipar )
323
       call matvecsubr( T1V, R, ipar )
324
       T1V = R - B
325
       call pcondlsubr( R, T1V, ipar )
326
       residual = normfun( HUTI_NDIM, R, 1 )
327
    end select
328

329
    work(m+1,S_ind) = normfun( HUTI_NDIM, R, 1 )
330

331
    !
332
    ! Print debugging output if desired
333
    !
334

335
    if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
336
       if ( mod(iter_count, HUTI_DBUGLVL) .eq. 0 ) then
337
          WRITE (*, '(A, I8, E11.4)') '   gmres:',iter_count, residual
338
       end if
339
    end if
340

341
    if ( residual .lt. HUTI_TOLERANCE ) then
342
       HUTI_INFO = HUTI_CONVERGENCE
343
       go to 1000
344
    end if
345

346
    IF( residual /= residual .OR. residual > HUTI_MAXTOLERANCE ) THEN
347
      HUTI_INFO = HUTI_DIVERGENCE
348
      GOTO 1000
349
    END IF
350

351
    !
352
    ! Return back to the iteration loop (without initialization)
353
    !
354

355
    iter_count = iter_count + 1
356
    if ( iter_count .gt. HUTI_MAXIT ) then
357
       HUTI_INFO = HUTI_MAXITER
358
       go to 1000
359
    end if
360

361
    go to 300
362

363
    !
364
    ! This is where we exit last time (after enough iterations or breakdown)
365
    !
366

367
1000 continue
368
    if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
369
       WRITE (*, '(A, I8, E11.4)') '   gmres:',iter_count, residual
370
    end if
371

372
    HUTI_ITERS = iter_count
373
    call pcondrsubr( T1V, X, ipar )
374
    X = T1V
375
    return
376

377
    ! End of execution
378
    !*********************************************************************
379

380
  end subroutine  huti_sgmressolv
381

382
  !*************************************************************************
383

384
  !*************************************************************************
385
  !*************************************************************************
386
  ! Double precision version
387
  !*************************************************************************
388
  !*************************************************************************
389

390
  subroutine  huti_dgmressolv  ( ndim, wrkdim, xvec, rhsvec, &
391
       ipar, dpar, work, matvecsubr, pcondlsubr, &
392
       pcondrsubr, dotprodfun, normfun, stopcfun )
393

394
    use huti_interfaces
395
    implicit none
396

397
    procedure( mv_iface_d ), pointer :: matvecsubr
398
    procedure( pc_iface_d ), pointer :: pcondlsubr
399
    procedure( pc_iface_d ), pointer :: pcondrsubr
400
    procedure( dotp_iface_d ), pointer :: dotprodfun
401
    procedure( norm_iface_d ), pointer :: normfun
402
    procedure( stopc_iface_d ), pointer :: stopcfun
403

404
    ! Parameters
405

406
    integer :: ndim, wrkdim
407
    double precision, dimension(ndim) :: xvec, rhsvec
408
    integer, dimension(HUTI_IPAR_DFLTSIZE) :: ipar
409
    double precision, dimension(HUTI_DPAR_DFLTSIZE) :: dpar
410
    double precision, dimension(ndim,wrkdim) :: work
411

412
    ! Local variables
413

414
    integer :: iter_count
415
    double precision :: residual, rhsnorm, precrhsnorm
416
    double precision :: bnrm, alpha, beta
417

418
    INTEGER :: reit, info, i, j, k, l, m, n, ii, jj
419

420
    double precision :: temp, temp2, error, gamma
421

422
    ! Local arrays
423

424
    double precision, DIMENSION(HUTI_GMRES_RESTART+1,HUTI_GMRES_RESTART+1) :: H
425
    double precision, &
426
         DIMENSION((HUTI_GMRES_RESTART+1)*(HUTI_GMRES_RESTART+1)) :: HLU
427
    double precision, DIMENSION(HUTI_GMRES_RESTART+1) :: CS, SN, Y
428

429
    !
430
    ! End of variable declarations
431
    !*********************************************************************
432

433
    !*********************************************************************
434
    ! The actual GMRES begins here (look the pseudo code in the
435
    ! "Templates..."-book, page 20)
436
    !
437
    ! First the initialization part
438
    !
439

440
    iter_count = 1
441

442
    bnrm = normfun( HUTI_NDIM, B, 1 )
443

444
    ! Norms of right-hand side vector are used in convergence tests
445

446
    if ( HUTI_STOPC .eq. HUTI_TRESID_SCALED_BYB .or. & 
447
         HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYB ) then
448
       rhsnorm = bnrm
449
    end if
450
    if ( HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYPRECB ) then
451
       call pcondlsubr( T1V, B, ipar )
452
       precrhsnorm = normfun( HUTI_NDIM, T1V, 1 )
453
    end if
454

455
    ! The following applies for all matrix operations in this solver
456

457
    HUTI_EXTOP_MATTYPE = HUTI_MAT_NOTTRPSED
458

459
    ! Generate vector X if needed
460

461
    if ( HUTI_INITIALX .eq. HUTI_RANDOMX ) then
462
       call  huti_drandvec   ( X, ipar )
463
    else if ( HUTI_INITIALX .ne. HUTI_USERSUPPLIEDX ) then
464
#ifdef _OPENMP
465
       !$OMP PARALLEL DO
466
       DO ii=1,ndim
467
         xvec(ii) = 1
468
       END DO
469
       !$OMP END PARALLEL DO
470
#else   
471
        X = 1
472
#endif
473
    end if
474

475
    call pcondrsubr( T1V, X, ipar )
476
    call matvecsubr( T1V, R, ipar )
477
#ifdef _OPENMP
478
    !$OMP PARALLEL DO
479
    DO ii=1,ndim
480
      work(ii,T1V_ind) = rhsvec(ii) - work(ii,R_ind)
481
    END DO
482
    !$OMP END PARALLEL DO
483
#else
484
    T1V = B - R
485
#endif
486
    call pcondlsubr( R, T1V, ipar )
487

488
    m = HUTI_GMRES_RESTART
489
#ifdef _OPENMP
490
    !$OMP PARALLEL PRIVATE(j)
491
    DO jj=V_ind+1-1,V_ind+m+1-1
492
      !$OMP DO
493
      DO ii=1,ndim
494
        work(ii,jj) = 0
495
      END DO
496
      !$OMP END DO
497
    END DO
498
    !$OMP DO
499
    DO ii=1,ndim
500
      work(ii,VTMP_ind) = 0
501
    END DO
502
    !$OMP END DO
503
    !$OMP END PARALLEL
504
#else
505
    work(:,V_ind+1-1:V_ind+m+1-1) = 0
506
    VTMP = 0
507
#endif
508
    H = 0
509
    CS = 0
510
    SN = 0
511
    work(1,VTMP_ind) = 1.0
512

513
    !
514
    ! This is where the loop starts (that is we continue from here after
515
    ! the first iteration)
516
    !
517

518
300 continue
519

520
    call pcondrsubr( T1V, X, ipar )
521
    call matvecsubr( T1V, R, ipar )
522
#ifdef _OPENMP
523
    !$OMP PARALLEL DO
524
    DO ii=1,ndim
525
      work(ii,T1V_ind) = rhsvec(ii) - work(ii,R_ind)
526
    END DO
527
    !$OMP END PARALLEL DO
528
#else
529
    T1V = B - R
530
#endif
531
    call pcondlsubr( R, T1V, ipar )
532

533
    alpha = normfun( HUTI_NDIM, R, 1 )
534
    if ( alpha .eq. 0 ) then
535
       HUTI_INFO = HUTI_GMRES_ALPHA
536
       go to 1000
537
    end if
538

539
#ifdef _OPENMP
540
    !$OMP PARALLEL
541
    !$OMP DO
542
    DO ii=1,ndim
543
      work(ii,V_ind+1-1) = work(ii,R_ind) / alpha
544
    END DO
545
    !$OMP END DO
546
    !$OMP DO
547
    DO ii=1,ndim
548
      work(ii,S_ind) = alpha * work(ii,VTMP_ind)
549
    END DO
550
    !$OMP END DO
551
    !$OMP END PARALLEL
552
#else
553
    work(:,V_ind+1-1) = R / alpha
554
    S = alpha * VTMP
555
#endif
556

557
    !
558
    ! Construct orthonormal
559
    !
560

561
    DO i = 1, m
562
       call pcondrsubr( W, work(:,V_ind+i-1), ipar )
563
       call matvecsubr( W, T1V, ipar )
564
       call pcondlsubr( W, T1V, ipar )
565

566
       DO k = 1, i
567
          H(k,i) = dotprodfun( HUTI_NDIM, W, 1, work(:,V_ind+k-1), 1 )
568
#ifdef _OPENMP
569
          !$OMP PARALLEL DO
570
          DO ii=1,ndim
571
            work(ii,W_ind) = work(ii,W_ind) - H(k,i) * work(ii,V_ind+k-1)
572
          END DO
573
          !$OMP END PARALLEL DO
574
#else
575
          W = W - H(k,i) * work(:,V_ind+k-1)
576
#endif
577
       END DO
578

579
       beta = normfun( HUTI_NDIM, W, 1 )
580
       if ( beta .eq. 0 ) then
581
          HUTI_INFO = HUTI_GMRES_BETA
582
          go to 1000
583
       end if
584

585
       H(i+1,i) = beta
586
#if _OPENMP
587
       !$OMP PARALLEL DO
588
       DO ii=1,ndim
589
         work(ii,V_ind+i+1-1) = work(ii,W_ind) / H(i+1, i)
590
       END DO
591
       !$OMP END PARALLEL DO
592
#else
593
       work(:,V_ind+i+1-1) = W / H(i+1, i)
594
#endif
595
       
596
       !
597
       ! Compute the Givens rotation
598
       !
599

600
       DO k = 1, i-1
601
          temp = CS(k) * H(k,i) + SN(k) * H(k+1,i)
602
          H(k+1,i) = -1 * SN(k) * H(k,i) + CS(k) * H(k+1,i)
603
          H(k,i) = temp
604
       END DO
605

606
       IF ( H(i+1,i) .eq. 0 ) THEN
607
          CS(i) = 1; SN(i) = 0
608
       ELSE
609
          IF ( abs( H(i+1,i) ) .gt. abs( H(i,i) ) ) THEN
610
             temp2 = H(i,i) / H(i+1,i)
611
             SN(i) = 1 / sqrt( 1 + ( temp2 * temp2 ))
612
             CS(i) = temp2 * SN(i)
613
          ELSE
614
             temp2 = H(i+1,i) / H(i,i)
615
             CS(i) = 1 / sqrt( 1 + ( temp2 * temp2 ))
616
             SN(i) = temp2 * CS(i)
617
          END IF
618
       END IF
619

620
       temp = CS(i) * work(i,S_ind)
621
       work(i+1,S_ind) = -1 * SN(i) * work(i,S_ind)
622
       work(i,S_ind) = temp
623

624
       H(i,i) = ( CS(i) * H(i,i) ) + ( SN(i) * H(i+1,i) )
625
       H(i+1,i) = 0
626

627
       error = abs( work(i+1,S_ind) ) / bnrm 
628
       IF ( REAL( error ) .lt.  HUTI_TOLERANCE ) THEN
629

630
          HLU = 0; j = 1
631
          do k = 1, i
632
             do l = 1, i
633
                HLU(j) = H(l,k)
634
                j = j + 1
635
             end do
636
          end do
637

638
          call  huti_dlusolve  ( i, HLU, Y, work(:,S_ind) )
639

640
#ifdef _OPENMP
641
          CALL DGEMV('N',ndim,i,1D0,work(1,V_ind),ndim,Y,1,1D0,xvec,1)
642
#else
643
          X = X + MATMUL( work(:,V_ind+1-1:V_ind+i-1), Y(1:i) )
644
#endif
645
          
646
          EXIT
647
       END IF
648

649
    END DO
650

651
    IF ( REAL( error ) .lt. HUTI_TOLERANCE ) THEN
652
       GOTO 500
653
    END IF
654

655
    HLU = 0; j = 1
656
    do k = 1, m
657
       do l = 1, m
658
          HLU(j) = H(l,k)
659
          j = j + 1
660
       end do
661
    end do
662

663
    call  huti_dlusolve  ( m, HLU, Y, work(:,S_ind) )
664

665
#ifdef _OPENMP
666
    CALL DGEMV('N',ndim,m,1D0,work(1,V_ind),ndim,Y,1,1D0,xvec,1)
667
#else
668
    X = X + MATMUL( work(:,V_ind+1-1:V_ind+m-1), Y(1:m) )
669
#endif
670

671

672
500 CONTINUE
673

674
    !
675
    ! Check the convergence
676
    !
677

678
    select case (HUTI_STOPC)
679
    case (HUTI_TRUERESIDUAL)
680
       call pcondrsubr( T1V, X, ipar )
681
       call matvecsubr( T1V, R, ipar )
682
#ifdef _OPENMP
683
       !$OMP PARALLEL DO
684
       DO ii=1,ndim
685
         work(ii,T1V_ind) = rhsvec(ii) - work(ii,R_ind)
686
       END DO
687
       !$OMP END PARALLEL DO
688
#else
689
       T1V = B - R
690
#endif
691
       call pcondlsubr( R, T1V, ipar )
692
       residual = normfun( HUTI_NDIM, R, 1 )
693
    case (HUTI_TRESID_SCALED_BYB)
694
       call pcondrsubr( T1V, X, ipar )
695
       call matvecsubr( T1V, R, ipar )
696
#ifdef _OPENMP
697
       !$OMP PARALLEL DO
698
       DO ii=1,ndim
699
         work(ii,T1V_ind) = rhsvec(ii) - work(ii,R_ind)
700
       END DO
701
       !$OMP END PARALLEL DO
702
#else
703
       T1V = B - R
704
#endif
705
       call pcondlsubr( R, T1V, ipar )
706
       residual = normfun( HUTI_NDIM, R, 1 ) / rhsnorm
707
    case (HUTI_PSEUDORESIDUAL)
708
       call matvecsubr( X, R, ipar )
709
#ifdef _OPENMP
710
       !$OMP PARALLEL DO
711
       DO ii=1,ndim
712
         work(ii,R_ind) = work(ii,R_ind) - rhsvec(ii)
713
       END DO
714
       !$OMP END PARALLEL DO
715
#else
716
       R = R - B
717
#endif
718
       call pcondlsubr( T1V, R, ipar )
719
       residual = normfun( HUTI_NDIM, T1V, 1 )
720
    case (HUTI_PRESID_SCALED_BYB)
721
       call matvecsubr( X, R, ipar )
722
#ifdef _OPENMP
723
       !$OMP PARALLEL DO
724
       DO ii=1,ndim
725
         work(ii,R_ind) = work(ii,R_ind) - rhsvec(ii)
726
       END DO
727
       !$OMP END PARALLEL DO
728
#else
729
       R = R - B
730
#endif
731
       call pcondlsubr( T1V, R, ipar )
732
       residual = normfun( HUTI_NDIM, T1V, 1 ) / rhsnorm
733
    case (HUTI_PRESID_SCALED_BYPRECB)
734
       call matvecsubr( X, R, ipar )
735
#ifdef _OPENMP
736
       !$OMP PARALLEL DO
737
       DO ii=1,ndim
738
         work(ii,R_ind) = work(ii,R_ind) - rhsvec(ii)
739
       END DO
740
       !$OMP END PARALLEL DO
741
#else
742
       R = R - B
743
#endif
744
       call pcondlsubr( T1V, R, ipar )
745
       residual = normfun( HUTI_NDIM, T1V, 1 ) / precrhsnorm
746
    case (HUTI_XDIFF_NORM)
747
#ifdef _OPENMP
748
       CALL DGEMV('N',ndim,m,1D0,work(1,V_ind),ndim,Y,1,0D0,work(1,T1V_ind),1)
749
#else
750
       T1V = MATMUL( work(:,V_ind+1-1:V_ind+m-1), Y(1:m) )
751
#endif
752
       residual = normfun( HUTI_NDIM, T1V, 1 )
753
    case (HUTI_USUPPLIED_STOPC)
754
       residual = stopcfun( X, B, R, ipar, dpar )
755
    case default
756
       call pcondrsubr( T1V, X, ipar )
757
       call matvecsubr( T1V, R, ipar )
758
#ifdef _OPENMP
759
       !$OMP PARALLEL DO
760
       DO ii=1,ndim
761
         work(ii,T1V_ind) = work(ii,R_ind) - rhsvec(ii)
762
       END DO
763
       !$OMP END PARALLEL DO
764
#else
765
       T1V = R - B
766
#endif
767
       call pcondlsubr( R, T1V, ipar )
768
       residual = normfun( HUTI_NDIM, R, 1 )
769
    end select
770

771
    work(m+1,S_ind) = normfun( HUTI_NDIM, R, 1 )
772

773
    !
774
    ! Print debugging output if desired
775
    !
776

777
    if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
778
       if ( mod(iter_count, HUTI_DBUGLVL) .eq. 0 ) then
779
          WRITE (*, '(A, I8, E11.4)') '   gmres:',iter_count, residual
780
       end if
781
    end if
782

783
    if ( residual .lt. HUTI_TOLERANCE ) then
784
       HUTI_INFO = HUTI_CONVERGENCE
785
       go to 1000
786
    end if
787

788
    IF( residual /= residual .OR. residual > HUTI_MAXTOLERANCE ) THEN
789
      HUTI_INFO = HUTI_DIVERGENCE
790
      GOTO 1000
791
    END IF
792

793
    !
794
    ! Return back to the iteration loop (without initialization)
795
    !
796

797
    iter_count = iter_count + 1
798
    if ( iter_count .gt. HUTI_MAXIT ) then
799
       HUTI_INFO = HUTI_MAXITER
800
       go to 1000
801
    end if
802

803
    go to 300
804

805
    !
806
    ! This is where we exit last time (after enough iterations or breakdown)
807
    !
808

809
1000 continue
810
    if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
811
       WRITE (*, '(A, I8, E11.4)') '   gmres:',iter_count, residual
812
    end if
813

814
    HUTI_ITERS = iter_count
815
    call pcondrsubr( T1V, X, ipar )
816
    X = T1V
817
    return
818

819
    ! End of execution
820
    !*********************************************************************
821

822
  end subroutine  huti_dgmressolv
823

824
  !*************************************************************************
825

826

827
  !*************************************************************************
828
  !*************************************************************************
829
  ! Complex version
830
  !*************************************************************************
831
  !*************************************************************************
832

833
  subroutine  huti_cgmressolv  ( ndim, wrkdim, xvec, rhsvec, &
834
       ipar, dpar, work, matvecsubr, pcondlsubr, &
835
       pcondrsubr, dotprodfun, normfun, stopcfun )
836

837
    use huti_interfaces
838
    implicit none
839

840
    procedure( mv_iface_c ), pointer :: matvecsubr
841
    procedure( pc_iface_c ), pointer :: pcondlsubr
842
    procedure( pc_iface_c ), pointer :: pcondrsubr
843
    procedure( dotp_iface_c ), pointer :: dotprodfun
844
    procedure( norm_iface_c ), pointer :: normfun
845
    procedure( stopc_iface_c ), pointer :: stopcfun
846

847
    ! Parameters
848

849
    integer :: ndim, wrkdim
850
    complex, dimension(ndim) :: xvec, rhsvec
851
    integer, dimension(HUTI_IPAR_DFLTSIZE) :: ipar
852
    double precision, dimension(HUTI_DPAR_DFLTSIZE) :: dpar
853
    complex, dimension(ndim,wrkdim) :: work
854

855
    ! Local variables
856

857
    integer :: iter_count
858
    real :: residual, rhsnorm, precrhsnorm
859
    real :: bnrm, alpha, beta
860

861
    INTEGER :: reit, info, i, j, k, l, m, n
862

863
    complex :: temp, temp2, error, gamma
864

865
    ! Local arrays
866

867
    complex, DIMENSION(HUTI_GMRES_RESTART+1,HUTI_GMRES_RESTART+1) :: H
868
    complex, &
869
         DIMENSION((HUTI_GMRES_RESTART+1)*(HUTI_GMRES_RESTART+1)) :: HLU
870
    complex, DIMENSION(HUTI_GMRES_RESTART+1) :: CS, SN, Y
871

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

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

883
    iter_count = 1
884

885
    bnrm = normfun( HUTI_NDIM, B, 1 )
886

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

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

898
    ! The following applies for all matrix operations in this solver
899

900
    HUTI_EXTOP_MATTYPE = HUTI_MAT_NOTTRPSED
901

902
    ! Generate vector X if needed
903

904
    if ( HUTI_INITIALX .eq. HUTI_RANDOMX ) then
905
       call  huti_crandvec   ( X, ipar )
906
    else if ( HUTI_INITIALX .ne. HUTI_USERSUPPLIEDX ) then
907
       X = 1
908
    end if
909

910
    call pcondrsubr( T1V, X, ipar )
911
    call matvecsubr( T1V, R, ipar )
912
    T1V = B - R
913
    call pcondlsubr( R, T1V, ipar )
914

915
    m = HUTI_GMRES_RESTART
916
    work(:,V_ind+1-1:V_ind+m+1-1) = 0
917
    H = 0
918
    CS = 0
919
    SN = 0
920
    VTMP = 0
921
    work(1,VTMP_ind) = 1.0
922

923
    !
924
    ! This is where the loop starts (that is we continue from here after
925
    ! the first iteration)
926
    !
927

928
300 continue
929

930
    call pcondrsubr( T1V, X, ipar )
931
    call matvecsubr( T1V, R, ipar )
932
    T1V = B - R
933
    call pcondlsubr( R, T1V, ipar )
934

935
    alpha = normfun( HUTI_NDIM, R, 1 )
936
    if ( alpha .eq. 0 ) then
937
       HUTI_INFO = HUTI_GMRES_ALPHA
938
       go to 1000
939
    end if
940

941
    work(:,V_ind+1-1) = R / alpha
942
    S = alpha * VTMP
943

944
    !
945
    ! Construct orthonormal
946
    !
947

948
    DO i = 1, m
949
       call pcondrsubr( W, work(:,V_ind+i-1), ipar )
950
       call matvecsubr( W, T1V, ipar )
951
       call pcondlsubr( W, T1V, ipar )
952

953
       DO k = 1, i
954
          H(k,i) = dotprodfun( HUTI_NDIM, W, 1, work(:,V_ind+k-1), 1 )
955
          W = W - H(k,i) * work(:,V_ind+k-1)
956
       END DO
957

958
       beta = normfun( HUTI_NDIM, W, 1 )
959
       if ( beta .eq. 0 ) then
960
          HUTI_INFO = HUTI_GMRES_BETA
961
          go to 1000
962
       end if
963

964
       H(i+1,i) = beta
965
       work(:,V_ind+i+1-1) = W / H(i+1, i)
966

967
       !
968
       ! Compute the Givens rotation
969
       !
970

971
       DO k = 1, i-1
972
          temp = CS(k) * H(k,i) + SN(k) * H(k+1,i)
973
          H(k+1,i) = -1 * SN(k) * H(k,i) + CS(k) * H(k+1,i)
974
          H(k,i) = temp
975
       END DO
976

977
       IF ( H(i+1,i) .eq. 0 ) THEN
978
          CS(i) = 1; SN(i) = 0
979
       ELSE
980
          IF ( abs( H(i+1,i) ) .gt. abs( H(i,i) ) ) THEN
981
             temp2 = H(i,i) / H(i+1,i)
982
             SN(i) = 1 / sqrt( 1 + ( temp2 * temp2 ))
983
             CS(i) = temp2 * SN(i)
984
          ELSE
985
             temp2 = H(i+1,i) / H(i,i)
986
             CS(i) = 1 / sqrt( 1 + ( temp2 * temp2 ))
987
             SN(i) = temp2 * CS(i)
988
          END IF
989
       END IF
990

991
       temp = CS(i) * work(i,S_ind)
992
       work(i+1,S_ind) = -1 * SN(i) * work(i,S_ind)
993
       work(i,S_ind) = temp
994

995
       H(i,i) = ( CS(i) * H(i,i) ) + ( SN(i) * H(i+1,i) )
996
       H(i+1,i) = 0
997

998
       error = abs( work(i+1,S_ind) ) / bnrm 
999
       IF ( REAL( error ) .lt.  HUTI_TOLERANCE ) THEN
1000

1001
          HLU = 0; j = 1
1002
          do k = 1, i
1003
             do l = 1, i
1004
                HLU(j) = H(l,k)
1005
                j = j + 1
1006
             end do
1007
          end do
1008

1009
          call  huti_clusolve  ( i, HLU, Y, work(:,S_ind) )
1010

1011
          X = X + MATMUL( work(:,V_ind+1-1:V_ind+i-1), Y(1:i) )
1012

1013
          EXIT
1014
       END IF
1015

1016
    END DO
1017

1018
    IF ( REAL( error ) .lt. HUTI_TOLERANCE ) THEN
1019
       GOTO 500
1020
    END IF
1021

1022
    HLU = 0; j = 1
1023
    do k = 1, m
1024
       do l = 1, m
1025
          HLU(j) = H(l,k)
1026
          j = j + 1
1027
       end do
1028
    end do
1029

1030
    call  huti_clusolve  ( m, HLU, Y, work(:,S_ind) )
1031

1032
    X = X + MATMUL( work(:,V_ind+1-1:V_ind+m-1), Y(1:m) )
1033

1034
500 CONTINUE
1035

1036
    !
1037
    ! Check the convergence
1038
    !
1039

1040
    select case (HUTI_STOPC)
1041
    case (HUTI_TRUERESIDUAL)
1042
       call pcondrsubr( T1V, X, ipar )
1043
       call matvecsubr( T1V, R, ipar )
1044
       T1V = B - R
1045
       call pcondlsubr( R, T1V, ipar )
1046
       residual = normfun( HUTI_NDIM, R, 1 )
1047
    case (HUTI_TRESID_SCALED_BYB)
1048
       call pcondrsubr( T1V, X, ipar )
1049
       call matvecsubr( T1V, R, ipar )
1050
       T1V = B - R
1051
       call pcondlsubr( R, T1V, ipar )
1052
       residual = normfun( HUTI_NDIM, R, 1 ) / rhsnorm
1053
    case (HUTI_PSEUDORESIDUAL)
1054
       call matvecsubr( X, R, ipar )
1055
       R = R - B
1056
       call pcondlsubr( T1V, R, ipar )
1057
       residual = normfun( HUTI_NDIM, T1V, 1 )
1058
    case (HUTI_PRESID_SCALED_BYB)
1059
       call matvecsubr( X, R, ipar )
1060
       R = R - B
1061
       call pcondlsubr( T1V, R, ipar )
1062
       residual = normfun( HUTI_NDIM, T1V, 1 ) / rhsnorm
1063
    case (HUTI_PRESID_SCALED_BYPRECB)
1064
       call matvecsubr( X, R, ipar )
1065
       R = R - B
1066
       call pcondlsubr( T1V, R, ipar )
1067
       residual = normfun( HUTI_NDIM, T1V, 1 ) / precrhsnorm
1068
    case (HUTI_XDIFF_NORM)
1069
       T1V = MATMUL( work(:,V_ind+1-1:V_ind+m-1), Y(1:m) )
1070
       residual = normfun( HUTI_NDIM, T1V, 1 )
1071
    case (HUTI_USUPPLIED_STOPC)
1072
       residual = stopcfun( X, B, R, ipar, dpar )
1073
    case default
1074
       call pcondrsubr( T1V, X, ipar )
1075
       call matvecsubr( T1V, R, ipar )
1076
       T1V = R - B
1077
       call pcondlsubr( R, T1V, ipar )
1078
       residual = normfun( HUTI_NDIM, R, 1 )
1079
    end select
1080

1081
    work(m+1,S_ind) = normfun( HUTI_NDIM, R, 1 )
1082

1083
    !
1084
    ! Print debugging output if desired
1085
    !
1086

1087
    if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
1088
       if ( mod(iter_count, HUTI_DBUGLVL) .eq. 0 ) then
1089
          WRITE (*, '(A, I8, E11.4)') '   gmresz:',iter_count, residual
1090
       end if
1091
    end if
1092

1093
    if ( residual .lt. HUTI_TOLERANCE ) then
1094
       HUTI_INFO = HUTI_CONVERGENCE
1095
       go to 1000
1096
    end if
1097

1098
    IF( residual /= residual .OR. residual > HUTI_MAXTOLERANCE ) THEN
1099
      HUTI_INFO = HUTI_DIVERGENCE
1100
      GOTO 1000
1101
    END IF
1102
    
1103
    !
1104
    ! Return back to the iteration loop (without initialization)
1105
    !
1106

1107
    iter_count = iter_count + 1
1108
    if ( iter_count .gt. HUTI_MAXIT ) then
1109
       HUTI_INFO = HUTI_MAXITER
1110
       go to 1000
1111
    end if
1112

1113
    go to 300
1114

1115
    !
1116
    ! This is where we exit last time (after enough iterations or breakdown)
1117
    !
1118

1119
1000 continue
1120
    if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
1121
       WRITE (*, '(A, I8, E11.4)') '   gmresz:', iter_count, residual
1122
    end if
1123

1124
    HUTI_ITERS = iter_count
1125
    call pcondrsubr( T1V, X, ipar )
1126
    X = T1V
1127
    return
1128

1129
    ! End of execution
1130
    !*********************************************************************
1131

1132
  end subroutine  huti_cgmressolv
1133

1134
  !*************************************************************************
1135

1136

1137
  !*************************************************************************
1138
  !*************************************************************************
1139
  ! Double complex version
1140
  !*************************************************************************
1141
  !*************************************************************************
1142

1143
  subroutine  huti_zgmressolv  ( ndim, wrkdim, xvec, rhsvec, &
1144
       ipar, dpar, work, matvecsubr, pcondlsubr, &
1145
       pcondrsubr, dotprodfun, normfun, stopcfun )
1146

1147
    use huti_interfaces
1148
    implicit none
1149

1150
    procedure( mv_iface_z ), pointer :: matvecsubr
1151
    procedure( pc_iface_z ), pointer :: pcondlsubr
1152
    procedure( pc_iface_z ), pointer :: pcondrsubr
1153
    procedure( dotp_iface_z ), pointer :: dotprodfun
1154
    procedure( norm_iface_z ), pointer :: normfun
1155
    procedure( stopc_iface_z ), pointer :: stopcfun
1156

1157
    ! Parameters
1158

1159
    integer :: ndim, wrkdim
1160
    double complex, dimension(ndim) :: xvec, rhsvec
1161
    integer, dimension(HUTI_IPAR_DFLTSIZE) :: ipar
1162
    double precision, dimension(HUTI_DPAR_DFLTSIZE) :: dpar
1163
    double complex, dimension(ndim,wrkdim) :: work
1164

1165
    ! Local variables
1166

1167
    integer :: iter_count
1168
    double precision :: residual, rhsnorm, precrhsnorm
1169
    double precision :: bnrm, alpha, beta
1170

1171
    INTEGER :: reit, info, i, j, k, l, m, n
1172

1173
    double complex :: temp, temp2, error, gamma
1174

1175
    ! Local arrays
1176

1177
    double complex, DIMENSION(HUTI_GMRES_RESTART+1,HUTI_GMRES_RESTART+1) :: H
1178
    double complex, &
1179
         DIMENSION((HUTI_GMRES_RESTART+1)*(HUTI_GMRES_RESTART+1)) :: HLU
1180
    double complex, DIMENSION(HUTI_GMRES_RESTART+1) :: CS, SN, Y
1181

1182
    !
1183
    ! End of variable declarations
1184
    !*********************************************************************
1185

1186
    !*********************************************************************
1187
    ! The actual GMRES begins here (look the pseudo code in the
1188
    ! "Templates..."-book, page 20)
1189
    !
1190
    ! First the initialization part
1191
    !
1192

1193
    iter_count = 1
1194

1195
    bnrm = normfun( HUTI_NDIM, B, 1 )
1196

1197
    ! Norms of right-hand side vector are used in convergence tests
1198

1199
    if ( HUTI_STOPC .eq. HUTI_TRESID_SCALED_BYB .or. & 
1200
         HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYB ) then
1201
       rhsnorm = bnrm
1202
    end if
1203
    if ( HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYPRECB ) then
1204
       call pcondlsubr( T1V, B, ipar )
1205
       precrhsnorm = normfun( HUTI_NDIM, T1V, 1 )
1206
    end if
1207

1208
    ! The following applies for all matrix operations in this solver
1209

1210
    HUTI_EXTOP_MATTYPE = HUTI_MAT_NOTTRPSED
1211

1212
    ! Generate vector X if needed
1213

1214
    if ( HUTI_INITIALX .eq. HUTI_RANDOMX ) then
1215
       call  huti_zrandvec   ( X, ipar )
1216
    else if ( HUTI_INITIALX .ne. HUTI_USERSUPPLIEDX ) then
1217
       X = 1
1218
    end if
1219

1220
    call pcondrsubr( T1V, X, ipar )
1221
    call matvecsubr( T1V, R, ipar )
1222
    T1V = B - R
1223
    call pcondlsubr( R, T1V, ipar )
1224

1225
    m = HUTI_GMRES_RESTART
1226
    work(:,V_ind+1-1:V_ind+m+1-1) = 0
1227
    H = 0
1228
    CS = 0
1229
    SN = 0
1230
    VTMP = 0
1231
    work(1,VTMP_ind) = 1.0
1232

1233
    !
1234
    ! This is where the loop starts (that is we continue from here after
1235
    ! the first iteration)
1236
    !
1237

1238
300 continue
1239

1240
    call pcondrsubr( T1V, X, ipar )
1241
    call matvecsubr( T1V, R, ipar )
1242
    T1V = B - R
1243
    call pcondlsubr( R, T1V, ipar )
1244

1245
    alpha = normfun( HUTI_NDIM, R, 1 )
1246
    if ( alpha .eq. 0 ) then
1247
       HUTI_INFO = HUTI_GMRES_ALPHA
1248
       go to 1000
1249
    end if
1250

1251
    work(:,V_ind+1-1) = R / alpha
1252
    S = alpha * VTMP
1253

1254
    !
1255
    ! Construct orthonormal
1256
    !
1257

1258
    DO i = 1, m
1259
       call pcondrsubr( W, work(:,V_ind+i-1), ipar )
1260
       call matvecsubr( W, T1V, ipar )
1261
       call pcondlsubr( W, T1V, ipar )
1262

1263
       DO k = 1, i
1264
          H(k,i) = dotprodfun( HUTI_NDIM, W, 1, work(:,V_ind+k-1), 1 )
1265
          W = W - H(k,i) * work(:,V_ind+k-1)
1266
       END DO
1267

1268
       beta = normfun( HUTI_NDIM, W, 1 )
1269
       if ( beta .eq. 0 ) then
1270
          HUTI_INFO = HUTI_GMRES_BETA
1271
          go to 1000
1272
       end if
1273

1274
       H(i+1,i) = beta
1275
       work(:,V_ind+i+1-1) = W / H(i+1, i)
1276

1277
       !
1278
       ! Compute the Givens rotation
1279
       !
1280

1281
       DO k = 1, i-1
1282
          temp = CS(k) * H(k,i) + SN(k) * H(k+1,i)
1283
          H(k+1,i) = -1 * SN(k) * H(k,i) + CS(k) * H(k+1,i)
1284
          H(k,i) = temp
1285
       END DO
1286

1287
       IF ( H(i+1,i) .eq. 0 ) THEN
1288
          CS(i) = 1; SN(i) = 0
1289
       ELSE
1290
          IF ( abs( H(i+1,i) ) .gt. abs( H(i,i) ) ) THEN
1291
             temp2 = H(i,i) / H(i+1,i)
1292
             SN(i) = 1 / sqrt( 1 + ( temp2 * temp2 ))
1293
             CS(i) = temp2 * SN(i)
1294
          ELSE
1295
             temp2 = H(i+1,i) / H(i,i)
1296
             CS(i) = 1 / sqrt( 1 + ( temp2 * temp2 ))
1297
             SN(i) = temp2 * CS(i)
1298
          END IF
1299
       END IF
1300

1301
       temp = CS(i) * work(i,S_ind)
1302
       work(i+1,S_ind) = -1 * SN(i) * work(i,S_ind)
1303
       work(i,S_ind) = temp
1304

1305
       H(i,i) = ( CS(i) * H(i,i) ) + ( SN(i) * H(i+1,i) )
1306
       H(i+1,i) = 0
1307

1308
       error = abs( work(i+1,S_ind) ) / bnrm 
1309
       IF ( REAL( error ) .lt.  HUTI_TOLERANCE ) THEN
1310

1311
          HLU = 0; j = 1
1312
          do k = 1, i
1313
             do l = 1, i
1314
                HLU(j) = H(l,k)
1315
                j = j + 1
1316
             end do
1317
          end do
1318

1319
          call  huti_zlusolve  ( i, HLU, Y, work(:,S_ind) )
1320

1321
          X = X + MATMUL( work(:,V_ind+1-1:V_ind+i-1), Y(1:i) )
1322

1323
          EXIT
1324
       END IF
1325

1326
    END DO
1327

1328
    IF ( REAL( error ) .lt. HUTI_TOLERANCE ) THEN
1329
       GOTO 500
1330
    END IF
1331

1332
    HLU = 0; j = 1
1333
    do k = 1, m
1334
       do l = 1, m
1335
          HLU(j) = H(l,k)
1336
          j = j + 1
1337
       end do
1338
    end do
1339

1340
    call  huti_zlusolve  ( m, HLU, Y, work(:,S_ind) )
1341

1342
    X = X + MATMUL( work(:,V_ind+1-1:V_ind+m-1), Y(1:m) )
1343

1344
500 CONTINUE
1345

1346
    !
1347
    ! Check the convergence
1348
    !
1349

1350
    select case (HUTI_STOPC)
1351
    case (HUTI_TRUERESIDUAL)
1352
       call pcondrsubr( T1V, X, ipar )
1353
       call matvecsubr( T1V, R, ipar )
1354
       T1V = B - R
1355
       call pcondlsubr( R, T1V, ipar )
1356
       residual = normfun( HUTI_NDIM, R, 1 )
1357
    case (HUTI_TRESID_SCALED_BYB)
1358
       call pcondrsubr( T1V, X, ipar )
1359
       call matvecsubr( T1V, R, ipar )
1360
       T1V = B - R
1361
       call pcondlsubr( R, T1V, ipar )
1362
       residual = normfun( HUTI_NDIM, R, 1 ) / rhsnorm
1363
    case (HUTI_PSEUDORESIDUAL)
1364
       call matvecsubr( X, R, ipar )
1365
       R = R - B
1366
       call pcondlsubr( T1V, R, ipar )
1367
       residual = normfun( HUTI_NDIM, T1V, 1 )
1368
    case (HUTI_PRESID_SCALED_BYB)
1369
       call matvecsubr( X, R, ipar )
1370
       R = R - B
1371
       call pcondlsubr( T1V, R, ipar )
1372
       residual = normfun( HUTI_NDIM, T1V, 1 ) / rhsnorm
1373
    case (HUTI_PRESID_SCALED_BYPRECB)
1374
       call matvecsubr( X, R, ipar )
1375
       R = R - B
1376
       call pcondlsubr( T1V, R, ipar )
1377
       residual = normfun( HUTI_NDIM, T1V, 1 ) / precrhsnorm
1378
    case (HUTI_XDIFF_NORM)
1379
       T1V = MATMUL( work(:,V_ind+1-1:V_ind+m-1), Y(1:m) )
1380
       residual = normfun( HUTI_NDIM, T1V, 1 )
1381
    case (HUTI_USUPPLIED_STOPC)
1382
       residual = stopcfun( X, B, R, ipar, dpar )
1383
    case default
1384
       call pcondrsubr( T1V, X, ipar )
1385
       call matvecsubr( T1V, R, ipar )
1386
       T1V = R - B
1387
       call pcondlsubr( R, T1V, ipar )
1388
       residual = normfun( HUTI_NDIM, R, 1 )
1389
    end select
1390

1391
    work(m+1,S_ind) = normfun( HUTI_NDIM, R, 1 )
1392

1393
    !
1394
    ! Print debugging output if desired
1395
    !
1396

1397
    if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
1398
       if ( mod(iter_count, HUTI_DBUGLVL) .eq. 0 ) then
1399
          WRITE (*, '(A, I8, E11.4)') '   gmresz:',iter_count, residual
1400
       end if
1401
    end if
1402

1403
    if ( residual .lt. HUTI_TOLERANCE ) then
1404
       HUTI_INFO = HUTI_CONVERGENCE
1405
       go to 1000
1406
    end if
1407

1408
    IF( residual /= residual .OR. residual > HUTI_MAXTOLERANCE ) THEN
1409
      HUTI_INFO = HUTI_DIVERGENCE
1410
      GOTO 1000
1411
    END IF
1412

1413
    
1414
    !
1415
    ! Return back to the iteration loop (without initialization)
1416
    !
1417

1418
    iter_count = iter_count + 1
1419
    if ( iter_count .gt. HUTI_MAXIT ) then
1420
       HUTI_INFO = HUTI_MAXITER
1421
       go to 1000
1422
    end if
1423

1424
    go to 300
1425

1426
    !
1427
    ! This is where we exit last time (after enough iterations or breakdown)
1428
    !
1429

1430
1000 continue
1431
    if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
1432
       WRITE (*, '(A,I8, E11.4)') '   gmresz:',iter_count, residual
1433
    end if
1434

1435
    HUTI_ITERS = iter_count
1436
    call pcondrsubr( T1V, X, ipar )
1437
    X = T1V
1438
    return
1439

1440
    ! End of execution
1441
    !*********************************************************************
1442

1443
  end subroutine  huti_zgmressolv
1444

1445
  !*************************************************************************
1446

1447

1448
end module huti_gmres
1449

1450
Product

Resources

Company
