Revision 600e72da161b072eb4d755f3b9685a147d07695e authored by Peter Raback on 23 November 2023, 10:18:35 UTC, committed by Peter Raback on 23 November 2023, 10:18:35 UTC
1 parent d809081
huti_gmres.F90
module huti_gmres
use huti_aux
implicit none
!
! *
! * Elmer, A Finite Element Software for Multiphysical Problems
! *
! * Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland
! *
! * This library is free software; you can redistribute it and/or
! * modify it under the terms of the GNU Lesser General Public
! * License as published by the Free Software Foundation; either
! * version 2.1 of the License, or (at your option) any later version.
! *
! * This library is distributed in the hope that it will be useful,
! * but WITHOUT ANY WARRANTY; without even the implied warranty of
! * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
! * Lesser General Public License for more details.
! *
! * You should have received a copy of the GNU Lesser General Public
! * License along with this library (in file ../LGPL-2.1); if not, write
! * to the Free Software Foundation, Inc., 51 Franklin Street,
! * Fifth Floor, Boston, MA 02110-1301 USA
! *
! *****************************************************************************/
!
! Subroutines to implement Generalized Minimum Residual iterative method
!
#include "huti_fdefs.h"
!*************************************************************************
!*************************************************************************
!
! These subroutines are based on a book by Barret et al.:
! "Templates for the Solution of Linear Systems: Building Blocks for
! Iterative Methods", 1993.
!
! All matrix-vector operations are done externally, so we do not need
! to know about the matrix structure (sparse or dense). So has the
! memory allocation for the working arrays done also externally.
!*************************************************************************
! Work array is used in the following order:
! work(:,1) = z
! work(:,2) = p
! work(:,3) = q
! work(:,4) = r
!
! Definitions to make the code more understandable and to make it look
! like the pseudo code (these are common to all precisions)
!
#define X xvec
#define B rhsvec
#define W work(:,1)
#define W_ind 1
#define R work(:,2)
#define R_ind 2
#define S work(:,3)
#define S_ind 3
#define VTMP work(:,4)
#define VTMP_ind 4
#define T1V work(:,5)
#define T1V_ind 5
#define V work(:,6)
#define V_ind 6
contains
!*************************************************************************
!*************************************************************************
!*************************************************************************
! Single precision version
!*************************************************************************
!*************************************************************************
subroutine huti_sgmressolv ( ndim, wrkdim, xvec, rhsvec, &
ipar, dpar, work, matvecsubr, pcondlsubr, &
pcondrsubr, dotprodfun, normfun, stopcfun )
use huti_interfaces
implicit none
procedure( mv_iface_s ), pointer :: matvecsubr
procedure( pc_iface_s ), pointer :: pcondlsubr
procedure( pc_iface_s ), pointer :: pcondrsubr
procedure( dotp_iface_s ), pointer :: dotprodfun
procedure( norm_iface_s ), pointer :: normfun
procedure( stopc_iface_s ), pointer :: stopcfun
! Parameters
integer :: ndim, wrkdim
real, dimension(ndim) :: xvec, rhsvec
integer, dimension(HUTI_IPAR_DFLTSIZE) :: ipar
double precision, dimension(HUTI_DPAR_DFLTSIZE) :: dpar
real, dimension(ndim,wrkdim) :: work
! Local variables
integer :: iter_count
real :: residual, rhsnorm, precrhsnorm
real :: bnrm, alpha, beta
INTEGER :: reit, info, i, j, k, l, m, n
real :: temp, temp2, error, gamma
! Local arrays
real, DIMENSION(HUTI_GMRES_RESTART+1,HUTI_GMRES_RESTART+1) :: H
real, &
DIMENSION((HUTI_GMRES_RESTART+1)*(HUTI_GMRES_RESTART+1)) :: HLU
real, DIMENSION(HUTI_GMRES_RESTART+1) :: CS, SN, Y
!
! End of variable declarations
!*********************************************************************
!*********************************************************************
! The actual GMRES begins here (look the pseudo code in the
! "Templates..."-book, page 20)
!
! First the initialization part
!
iter_count = 1
bnrm = normfun( HUTI_NDIM, B, 1 )
! Norms of right-hand side vector are used in convergence tests
if ( HUTI_STOPC .eq. HUTI_TRESID_SCALED_BYB .or. &
HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYB ) then
rhsnorm = bnrm
end if
if ( HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYPRECB ) then
call pcondlsubr( T1V, B, ipar )
precrhsnorm = normfun( HUTI_NDIM, T1V, 1 )
end if
! The following applies for all matrix operations in this solver
HUTI_EXTOP_MATTYPE = HUTI_MAT_NOTTRPSED
! Generate vector X if needed
if ( HUTI_INITIALX .eq. HUTI_RANDOMX ) then
call huti_srandvec ( X, ipar )
else if ( HUTI_INITIALX .ne. HUTI_USERSUPPLIEDX ) then
X = 1
end if
call pcondrsubr( T1V, X, ipar )
call matvecsubr( T1V, R, ipar )
T1V = B - R
call pcondlsubr( R, T1V, ipar )
m = HUTI_GMRES_RESTART
work(:,V_ind+1-1:V_ind+m+1-1) = 0
H = 0
CS = 0
SN = 0
VTMP = 0
work(1,VTMP_ind) = 1.0
!
! This is where the loop starts (that is we continue from here after
! the first iteration)
!
300 continue
call pcondrsubr( T1V, X, ipar )
call matvecsubr( T1V, R, ipar )
T1V = B - R
call pcondlsubr( R, T1V, ipar )
alpha = normfun( HUTI_NDIM, R, 1 )
if ( alpha .eq. 0 ) then
HUTI_INFO = HUTI_GMRES_ALPHA
go to 1000
end if
work(:,V_ind+1-1) = R / alpha
S = alpha * VTMP
!
! Construct orthonormal
!
DO i = 1, m
call pcondrsubr( W, work(:,V_ind+i-1), ipar )
call matvecsubr( W, T1V, ipar )
call pcondlsubr( W, T1V, ipar )
DO k = 1, i
H(k,i) = dotprodfun( HUTI_NDIM, W, 1, work(:,V_ind+k-1), 1 )
W = W - H(k,i) * work(:,V_ind+k-1)
END DO
beta = normfun( HUTI_NDIM, W, 1 )
if ( beta .eq. 0 ) then
HUTI_INFO = HUTI_GMRES_BETA
go to 1000
end if
H(i+1,i) = beta
work(:,V_ind+i+1-1) = W / H(i+1, i)
!
! Compute the Givens rotation
!
DO k = 1, i-1
temp = CS(k) * H(k,i) + SN(k) * H(k+1,i)
H(k+1,i) = -1 * SN(k) * H(k,i) + CS(k) * H(k+1,i)
H(k,i) = temp
END DO
IF ( H(i+1,i) .eq. 0 ) THEN
CS(i) = 1; SN(i) = 0
ELSE
IF ( abs( H(i+1,i) ) .gt. abs( H(i,i) ) ) THEN
temp2 = H(i,i) / H(i+1,i)
SN(i) = 1 / sqrt( 1 + ( temp2 * temp2 ))
CS(i) = temp2 * SN(i)
ELSE
temp2 = H(i+1,i) / H(i,i)
CS(i) = 1 / sqrt( 1 + ( temp2 * temp2 ))
SN(i) = temp2 * CS(i)
END IF
END IF
temp = CS(i) * work(i,S_ind)
work(i+1,S_ind) = -1 * SN(i) * work(i,S_ind)
work(i,S_ind) = temp
H(i,i) = ( CS(i) * H(i,i) ) + ( SN(i) * H(i+1,i) )
H(i+1,i) = 0
error = abs( work(i+1,S_ind) ) / bnrm
IF ( REAL( error ) .lt. HUTI_TOLERANCE ) THEN
HLU = 0; j = 1
do k = 1, i
do l = 1, i
HLU(j) = H(l,k)
j = j + 1
end do
end do
call huti_slusolve ( i, HLU, Y, work(:,S_ind) )
X = X + MATMUL( work(:,V_ind+1-1:V_ind+i-1), Y(1:i) )
EXIT
END IF
END DO
IF ( REAL( error ) .lt. HUTI_TOLERANCE ) THEN
GOTO 500
END IF
HLU = 0; j = 1
do k = 1, m
do l = 1, m
HLU(j) = H(l,k)
j = j + 1
end do
end do
call huti_slusolve ( m, HLU, Y, work(:,S_ind) )
X = X + MATMUL( work(:,V_ind+1-1:V_ind+m-1), Y(1:m) )
500 CONTINUE
!
! Check the convergence
!
select case (HUTI_STOPC)
case (HUTI_TRUERESIDUAL)
call pcondrsubr( T1V, X, ipar )
call matvecsubr( T1V, R, ipar )
T1V = B - R
call pcondlsubr( R, T1V, ipar )
residual = normfun( HUTI_NDIM, R, 1 )
case (HUTI_TRESID_SCALED_BYB)
call pcondrsubr( T1V, X, ipar )
call matvecsubr( T1V, R, ipar )
T1V = B - R
call pcondlsubr( R, T1V, ipar )
residual = normfun( HUTI_NDIM, R, 1 ) / rhsnorm
case (HUTI_PSEUDORESIDUAL)
call matvecsubr( X, R, ipar )
R = R - B
call pcondlsubr( T1V, R, ipar )
residual = normfun( HUTI_NDIM, T1V, 1 )
case (HUTI_PRESID_SCALED_BYB)
call matvecsubr( X, R, ipar )
R = R - B
call pcondlsubr( T1V, R, ipar )
residual = normfun( HUTI_NDIM, T1V, 1 ) / rhsnorm
case (HUTI_PRESID_SCALED_BYPRECB)
call matvecsubr( X, R, ipar )
R = R - B
call pcondlsubr( T1V, R, ipar )
residual = normfun( HUTI_NDIM, T1V, 1 ) / precrhsnorm
case (HUTI_XDIFF_NORM)
T1V = MATMUL( work(:,V_ind+1-1:V_ind+m-1), Y(1:m) )
residual = normfun( HUTI_NDIM, T1V, 1 )
case (HUTI_USUPPLIED_STOPC)
residual = stopcfun( X, B, R, ipar, dpar )
case default
call pcondrsubr( T1V, X, ipar )
call matvecsubr( T1V, R, ipar )
T1V = R - B
call pcondlsubr( R, T1V, ipar )
residual = normfun( HUTI_NDIM, R, 1 )
end select
work(m+1,S_ind) = normfun( HUTI_NDIM, R, 1 )
!
! Print debugging output if desired
!
if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
if ( mod(iter_count, HUTI_DBUGLVL) .eq. 0 ) then
WRITE (*, '(A, I8, E11.4)') ' gmres:',iter_count, residual
end if
end if
if ( residual .lt. HUTI_TOLERANCE ) then
HUTI_INFO = HUTI_CONVERGENCE
go to 1000
end if
IF( residual /= residual .OR. residual > HUTI_MAXTOLERANCE ) THEN
HUTI_INFO = HUTI_DIVERGENCE
GOTO 1000
END IF
!
! Return back to the iteration loop (without initialization)
!
iter_count = iter_count + 1
if ( iter_count .gt. HUTI_MAXIT ) then
HUTI_INFO = HUTI_MAXITER
go to 1000
end if
go to 300
!
! This is where we exit last time (after enough iterations or breakdown)
!
1000 continue
if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
WRITE (*, '(A, I8, E11.4)') ' gmres:',iter_count, residual
end if
HUTI_ITERS = iter_count
call pcondrsubr( T1V, X, ipar )
X = T1V
return
! End of execution
!*********************************************************************
end subroutine huti_sgmressolv
!*************************************************************************
!*************************************************************************
!*************************************************************************
! Double precision version
!*************************************************************************
!*************************************************************************
subroutine huti_dgmressolv ( ndim, wrkdim, xvec, rhsvec, &
ipar, dpar, work, matvecsubr, pcondlsubr, &
pcondrsubr, dotprodfun, normfun, stopcfun )
use huti_interfaces
implicit none
procedure( mv_iface_d ), pointer :: matvecsubr
procedure( pc_iface_d ), pointer :: pcondlsubr
procedure( pc_iface_d ), pointer :: pcondrsubr
procedure( dotp_iface_d ), pointer :: dotprodfun
procedure( norm_iface_d ), pointer :: normfun
procedure( stopc_iface_d ), pointer :: stopcfun
! Parameters
integer :: ndim, wrkdim
double precision, dimension(ndim) :: xvec, rhsvec
integer, dimension(HUTI_IPAR_DFLTSIZE) :: ipar
double precision, dimension(HUTI_DPAR_DFLTSIZE) :: dpar
double precision, dimension(ndim,wrkdim) :: work
! Local variables
integer :: iter_count
double precision :: residual, rhsnorm, precrhsnorm
double precision :: bnrm, alpha, beta
INTEGER :: reit, info, i, j, k, l, m, n, ii, jj
double precision :: temp, temp2, error, gamma
! Local arrays
double precision, DIMENSION(HUTI_GMRES_RESTART+1,HUTI_GMRES_RESTART+1) :: H
double precision, &
DIMENSION((HUTI_GMRES_RESTART+1)*(HUTI_GMRES_RESTART+1)) :: HLU
double precision, DIMENSION(HUTI_GMRES_RESTART+1) :: CS, SN, Y
!
! End of variable declarations
!*********************************************************************
!*********************************************************************
! The actual GMRES begins here (look the pseudo code in the
! "Templates..."-book, page 20)
!
! First the initialization part
!
iter_count = 1
bnrm = normfun( HUTI_NDIM, B, 1 )
! Norms of right-hand side vector are used in convergence tests
if ( HUTI_STOPC .eq. HUTI_TRESID_SCALED_BYB .or. &
HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYB ) then
rhsnorm = bnrm
end if
if ( HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYPRECB ) then
call pcondlsubr( T1V, B, ipar )
precrhsnorm = normfun( HUTI_NDIM, T1V, 1 )
end if
! The following applies for all matrix operations in this solver
HUTI_EXTOP_MATTYPE = HUTI_MAT_NOTTRPSED
! Generate vector X if needed
if ( HUTI_INITIALX .eq. HUTI_RANDOMX ) then
call huti_drandvec ( X, ipar )
else if ( HUTI_INITIALX .ne. HUTI_USERSUPPLIEDX ) then
#ifdef _OPENMP
!$OMP PARALLEL DO
DO ii=1,ndim
xvec(ii) = 1
END DO
!$OMP END PARALLEL DO
#else
X = 1
#endif
end if
call pcondrsubr( T1V, X, ipar )
call matvecsubr( T1V, R, ipar )
#ifdef _OPENMP
!$OMP PARALLEL DO
DO ii=1,ndim
work(ii,T1V_ind) = rhsvec(ii) - work(ii,R_ind)
END DO
!$OMP END PARALLEL DO
#else
T1V = B - R
#endif
call pcondlsubr( R, T1V, ipar )
m = HUTI_GMRES_RESTART
#ifdef _OPENMP
!$OMP PARALLEL PRIVATE(j)
DO jj=V_ind+1-1,V_ind+m+1-1
!$OMP DO
DO ii=1,ndim
work(ii,jj) = 0
END DO
!$OMP END DO
END DO
!$OMP DO
DO ii=1,ndim
work(ii,VTMP_ind) = 0
END DO
!$OMP END DO
!$OMP END PARALLEL
#else
work(:,V_ind+1-1:V_ind+m+1-1) = 0
VTMP = 0
#endif
H = 0
CS = 0
SN = 0
work(1,VTMP_ind) = 1.0
!
! This is where the loop starts (that is we continue from here after
! the first iteration)
!
300 continue
call pcondrsubr( T1V, X, ipar )
call matvecsubr( T1V, R, ipar )
#ifdef _OPENMP
!$OMP PARALLEL DO
DO ii=1,ndim
work(ii,T1V_ind) = rhsvec(ii) - work(ii,R_ind)
END DO
!$OMP END PARALLEL DO
#else
T1V = B - R
#endif
call pcondlsubr( R, T1V, ipar )
alpha = normfun( HUTI_NDIM, R, 1 )
if ( alpha .eq. 0 ) then
HUTI_INFO = HUTI_GMRES_ALPHA
go to 1000
end if
#ifdef _OPENMP
!$OMP PARALLEL
!$OMP DO
DO ii=1,ndim
work(ii,V_ind+1-1) = work(ii,R_ind) / alpha
END DO
!$OMP END DO
!$OMP DO
DO ii=1,ndim
work(ii,S_ind) = alpha * work(ii,VTMP_ind)
END DO
!$OMP END DO
!$OMP END PARALLEL
#else
work(:,V_ind+1-1) = R / alpha
S = alpha * VTMP
#endif
!
! Construct orthonormal
!
DO i = 1, m
call pcondrsubr( W, work(:,V_ind+i-1), ipar )
call matvecsubr( W, T1V, ipar )
call pcondlsubr( W, T1V, ipar )
DO k = 1, i
H(k,i) = dotprodfun( HUTI_NDIM, W, 1, work(:,V_ind+k-1), 1 )
#ifdef _OPENMP
!$OMP PARALLEL DO
DO ii=1,ndim
work(ii,W_ind) = work(ii,W_ind) - H(k,i) * work(ii,V_ind+k-1)
END DO
!$OMP END PARALLEL DO
#else
W = W - H(k,i) * work(:,V_ind+k-1)
#endif
END DO
beta = normfun( HUTI_NDIM, W, 1 )
if ( beta .eq. 0 ) then
HUTI_INFO = HUTI_GMRES_BETA
go to 1000
end if
H(i+1,i) = beta
#if _OPENMP
!$OMP PARALLEL DO
DO ii=1,ndim
work(ii,V_ind+i+1-1) = work(ii,W_ind) / H(i+1, i)
END DO
!$OMP END PARALLEL DO
#else
work(:,V_ind+i+1-1) = W / H(i+1, i)
#endif
!
! Compute the Givens rotation
!
DO k = 1, i-1
temp = CS(k) * H(k,i) + SN(k) * H(k+1,i)
H(k+1,i) = -1 * SN(k) * H(k,i) + CS(k) * H(k+1,i)
H(k,i) = temp
END DO
IF ( H(i+1,i) .eq. 0 ) THEN
CS(i) = 1; SN(i) = 0
ELSE
IF ( abs( H(i+1,i) ) .gt. abs( H(i,i) ) ) THEN
temp2 = H(i,i) / H(i+1,i)
SN(i) = 1 / sqrt( 1 + ( temp2 * temp2 ))
CS(i) = temp2 * SN(i)
ELSE
temp2 = H(i+1,i) / H(i,i)
CS(i) = 1 / sqrt( 1 + ( temp2 * temp2 ))
SN(i) = temp2 * CS(i)
END IF
END IF
temp = CS(i) * work(i,S_ind)
work(i+1,S_ind) = -1 * SN(i) * work(i,S_ind)
work(i,S_ind) = temp
H(i,i) = ( CS(i) * H(i,i) ) + ( SN(i) * H(i+1,i) )
H(i+1,i) = 0
error = abs( work(i+1,S_ind) ) / bnrm
IF ( REAL( error ) .lt. HUTI_TOLERANCE ) THEN
HLU = 0; j = 1
do k = 1, i
do l = 1, i
HLU(j) = H(l,k)
j = j + 1
end do
end do
call huti_dlusolve ( i, HLU, Y, work(:,S_ind) )
#ifdef _OPENMP
CALL DGEMV('N',ndim,i,1D0,work(1,V_ind),ndim,Y,1,1D0,xvec,1)
#else
X = X + MATMUL( work(:,V_ind+1-1:V_ind+i-1), Y(1:i) )
#endif
EXIT
END IF
END DO
IF ( REAL( error ) .lt. HUTI_TOLERANCE ) THEN
GOTO 500
END IF
HLU = 0; j = 1
do k = 1, m
do l = 1, m
HLU(j) = H(l,k)
j = j + 1
end do
end do
call huti_dlusolve ( m, HLU, Y, work(:,S_ind) )
#ifdef _OPENMP
CALL DGEMV('N',ndim,m,1D0,work(1,V_ind),ndim,Y,1,1D0,xvec,1)
#else
X = X + MATMUL( work(:,V_ind+1-1:V_ind+m-1), Y(1:m) )
#endif
500 CONTINUE
!
! Check the convergence
!
select case (HUTI_STOPC)
case (HUTI_TRUERESIDUAL)
call pcondrsubr( T1V, X, ipar )
call matvecsubr( T1V, R, ipar )
#ifdef _OPENMP
!$OMP PARALLEL DO
DO ii=1,ndim
work(ii,T1V_ind) = rhsvec(ii) - work(ii,R_ind)
END DO
!$OMP END PARALLEL DO
#else
T1V = B - R
#endif
call pcondlsubr( R, T1V, ipar )
residual = normfun( HUTI_NDIM, R, 1 )
case (HUTI_TRESID_SCALED_BYB)
call pcondrsubr( T1V, X, ipar )
call matvecsubr( T1V, R, ipar )
#ifdef _OPENMP
!$OMP PARALLEL DO
DO ii=1,ndim
work(ii,T1V_ind) = rhsvec(ii) - work(ii,R_ind)
END DO
!$OMP END PARALLEL DO
#else
T1V = B - R
#endif
call pcondlsubr( R, T1V, ipar )
residual = normfun( HUTI_NDIM, R, 1 ) / rhsnorm
case (HUTI_PSEUDORESIDUAL)
call matvecsubr( X, R, ipar )
#ifdef _OPENMP
!$OMP PARALLEL DO
DO ii=1,ndim
work(ii,R_ind) = work(ii,R_ind) - rhsvec(ii)
END DO
!$OMP END PARALLEL DO
#else
R = R - B
#endif
call pcondlsubr( T1V, R, ipar )
residual = normfun( HUTI_NDIM, T1V, 1 )
case (HUTI_PRESID_SCALED_BYB)
call matvecsubr( X, R, ipar )
#ifdef _OPENMP
!$OMP PARALLEL DO
DO ii=1,ndim
work(ii,R_ind) = work(ii,R_ind) - rhsvec(ii)
END DO
!$OMP END PARALLEL DO
#else
R = R - B
#endif
call pcondlsubr( T1V, R, ipar )
residual = normfun( HUTI_NDIM, T1V, 1 ) / rhsnorm
case (HUTI_PRESID_SCALED_BYPRECB)
call matvecsubr( X, R, ipar )
#ifdef _OPENMP
!$OMP PARALLEL DO
DO ii=1,ndim
work(ii,R_ind) = work(ii,R_ind) - rhsvec(ii)
END DO
!$OMP END PARALLEL DO
#else
R = R - B
#endif
call pcondlsubr( T1V, R, ipar )
residual = normfun( HUTI_NDIM, T1V, 1 ) / precrhsnorm
case (HUTI_XDIFF_NORM)
#ifdef _OPENMP
CALL DGEMV('N',ndim,m,1D0,work(1,V_ind),ndim,Y,1,0D0,work(1,T1V_ind),1)
#else
T1V = MATMUL( work(:,V_ind+1-1:V_ind+m-1), Y(1:m) )
#endif
residual = normfun( HUTI_NDIM, T1V, 1 )
case (HUTI_USUPPLIED_STOPC)
residual = stopcfun( X, B, R, ipar, dpar )
case default
call pcondrsubr( T1V, X, ipar )
call matvecsubr( T1V, R, ipar )
#ifdef _OPENMP
!$OMP PARALLEL DO
DO ii=1,ndim
work(ii,T1V_ind) = work(ii,R_ind) - rhsvec(ii)
END DO
!$OMP END PARALLEL DO
#else
T1V = R - B
#endif
call pcondlsubr( R, T1V, ipar )
residual = normfun( HUTI_NDIM, R, 1 )
end select
work(m+1,S_ind) = normfun( HUTI_NDIM, R, 1 )
!
! Print debugging output if desired
!
if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
if ( mod(iter_count, HUTI_DBUGLVL) .eq. 0 ) then
WRITE (*, '(A, I8, E11.4)') ' gmres:',iter_count, residual
end if
end if
if ( residual .lt. HUTI_TOLERANCE ) then
HUTI_INFO = HUTI_CONVERGENCE
go to 1000
end if
IF( residual /= residual .OR. residual > HUTI_MAXTOLERANCE ) THEN
HUTI_INFO = HUTI_DIVERGENCE
GOTO 1000
END IF
!
! Return back to the iteration loop (without initialization)
!
iter_count = iter_count + 1
if ( iter_count .gt. HUTI_MAXIT ) then
HUTI_INFO = HUTI_MAXITER
go to 1000
end if
go to 300
!
! This is where we exit last time (after enough iterations or breakdown)
!
1000 continue
if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
WRITE (*, '(A, I8, E11.4)') ' gmres:',iter_count, residual
end if
HUTI_ITERS = iter_count
call pcondrsubr( T1V, X, ipar )
X = T1V
return
! End of execution
!*********************************************************************
end subroutine huti_dgmressolv
!*************************************************************************
!*************************************************************************
!*************************************************************************
! Complex version
!*************************************************************************
!*************************************************************************
subroutine huti_cgmressolv ( ndim, wrkdim, xvec, rhsvec, &
ipar, dpar, work, matvecsubr, pcondlsubr, &
pcondrsubr, dotprodfun, normfun, stopcfun )
use huti_interfaces
implicit none
procedure( mv_iface_c ), pointer :: matvecsubr
procedure( pc_iface_c ), pointer :: pcondlsubr
procedure( pc_iface_c ), pointer :: pcondrsubr
procedure( dotp_iface_c ), pointer :: dotprodfun
procedure( norm_iface_c ), pointer :: normfun
procedure( stopc_iface_c ), pointer :: stopcfun
! Parameters
integer :: ndim, wrkdim
complex, dimension(ndim) :: xvec, rhsvec
integer, dimension(HUTI_IPAR_DFLTSIZE) :: ipar
double precision, dimension(HUTI_DPAR_DFLTSIZE) :: dpar
complex, dimension(ndim,wrkdim) :: work
! Local variables
integer :: iter_count
real :: residual, rhsnorm, precrhsnorm
real :: bnrm, alpha, beta
INTEGER :: reit, info, i, j, k, l, m, n
complex :: temp, temp2, error, gamma
! Local arrays
complex, DIMENSION(HUTI_GMRES_RESTART+1,HUTI_GMRES_RESTART+1) :: H
complex, &
DIMENSION((HUTI_GMRES_RESTART+1)*(HUTI_GMRES_RESTART+1)) :: HLU
complex, DIMENSION(HUTI_GMRES_RESTART+1) :: CS, SN, Y
!
! End of variable declarations
!*********************************************************************
!*********************************************************************
! The actual GMRES begins here (look the pseudo code in the
! "Templates..."-book, page 20)
!
! First the initialization part
!
iter_count = 1
bnrm = normfun( HUTI_NDIM, B, 1 )
! Norms of right-hand side vector are used in convergence tests
if ( HUTI_STOPC .eq. HUTI_TRESID_SCALED_BYB .or. &
HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYB ) then
rhsnorm = bnrm
end if
if ( HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYPRECB ) then
call pcondlsubr( T1V, B, ipar )
precrhsnorm = normfun( HUTI_NDIM, T1V, 1 )
end if
! The following applies for all matrix operations in this solver
HUTI_EXTOP_MATTYPE = HUTI_MAT_NOTTRPSED
! Generate vector X if needed
if ( HUTI_INITIALX .eq. HUTI_RANDOMX ) then
call huti_crandvec ( X, ipar )
else if ( HUTI_INITIALX .ne. HUTI_USERSUPPLIEDX ) then
X = 1
end if
call pcondrsubr( T1V, X, ipar )
call matvecsubr( T1V, R, ipar )
T1V = B - R
call pcondlsubr( R, T1V, ipar )
m = HUTI_GMRES_RESTART
work(:,V_ind+1-1:V_ind+m+1-1) = 0
H = 0
CS = 0
SN = 0
VTMP = 0
work(1,VTMP_ind) = 1.0
!
! This is where the loop starts (that is we continue from here after
! the first iteration)
!
300 continue
call pcondrsubr( T1V, X, ipar )
call matvecsubr( T1V, R, ipar )
T1V = B - R
call pcondlsubr( R, T1V, ipar )
alpha = normfun( HUTI_NDIM, R, 1 )
if ( alpha .eq. 0 ) then
HUTI_INFO = HUTI_GMRES_ALPHA
go to 1000
end if
work(:,V_ind+1-1) = R / alpha
S = alpha * VTMP
!
! Construct orthonormal
!
DO i = 1, m
call pcondrsubr( W, work(:,V_ind+i-1), ipar )
call matvecsubr( W, T1V, ipar )
call pcondlsubr( W, T1V, ipar )
DO k = 1, i
H(k,i) = dotprodfun( HUTI_NDIM, W, 1, work(:,V_ind+k-1), 1 )
W = W - H(k,i) * work(:,V_ind+k-1)
END DO
beta = normfun( HUTI_NDIM, W, 1 )
if ( beta .eq. 0 ) then
HUTI_INFO = HUTI_GMRES_BETA
go to 1000
end if
H(i+1,i) = beta
work(:,V_ind+i+1-1) = W / H(i+1, i)
!
! Compute the Givens rotation
!
DO k = 1, i-1
temp = CS(k) * H(k,i) + SN(k) * H(k+1,i)
H(k+1,i) = -1 * SN(k) * H(k,i) + CS(k) * H(k+1,i)
H(k,i) = temp
END DO
IF ( H(i+1,i) .eq. 0 ) THEN
CS(i) = 1; SN(i) = 0
ELSE
IF ( abs( H(i+1,i) ) .gt. abs( H(i,i) ) ) THEN
temp2 = H(i,i) / H(i+1,i)
SN(i) = 1 / sqrt( 1 + ( temp2 * temp2 ))
CS(i) = temp2 * SN(i)
ELSE
temp2 = H(i+1,i) / H(i,i)
CS(i) = 1 / sqrt( 1 + ( temp2 * temp2 ))
SN(i) = temp2 * CS(i)
END IF
END IF
temp = CS(i) * work(i,S_ind)
work(i+1,S_ind) = -1 * SN(i) * work(i,S_ind)
work(i,S_ind) = temp
H(i,i) = ( CS(i) * H(i,i) ) + ( SN(i) * H(i+1,i) )
H(i+1,i) = 0
error = abs( work(i+1,S_ind) ) / bnrm
IF ( REAL( error ) .lt. HUTI_TOLERANCE ) THEN
HLU = 0; j = 1
do k = 1, i
do l = 1, i
HLU(j) = H(l,k)
j = j + 1
end do
end do
call huti_clusolve ( i, HLU, Y, work(:,S_ind) )
X = X + MATMUL( work(:,V_ind+1-1:V_ind+i-1), Y(1:i) )
EXIT
END IF
END DO
IF ( REAL( error ) .lt. HUTI_TOLERANCE ) THEN
GOTO 500
END IF
HLU = 0; j = 1
do k = 1, m
do l = 1, m
HLU(j) = H(l,k)
j = j + 1
end do
end do
call huti_clusolve ( m, HLU, Y, work(:,S_ind) )
X = X + MATMUL( work(:,V_ind+1-1:V_ind+m-1), Y(1:m) )
500 CONTINUE
!
! Check the convergence
!
select case (HUTI_STOPC)
case (HUTI_TRUERESIDUAL)
call pcondrsubr( T1V, X, ipar )
call matvecsubr( T1V, R, ipar )
T1V = B - R
call pcondlsubr( R, T1V, ipar )
residual = normfun( HUTI_NDIM, R, 1 )
case (HUTI_TRESID_SCALED_BYB)
call pcondrsubr( T1V, X, ipar )
call matvecsubr( T1V, R, ipar )
T1V = B - R
call pcondlsubr( R, T1V, ipar )
residual = normfun( HUTI_NDIM, R, 1 ) / rhsnorm
case (HUTI_PSEUDORESIDUAL)
call matvecsubr( X, R, ipar )
R = R - B
call pcondlsubr( T1V, R, ipar )
residual = normfun( HUTI_NDIM, T1V, 1 )
case (HUTI_PRESID_SCALED_BYB)
call matvecsubr( X, R, ipar )
R = R - B
call pcondlsubr( T1V, R, ipar )
residual = normfun( HUTI_NDIM, T1V, 1 ) / rhsnorm
case (HUTI_PRESID_SCALED_BYPRECB)
call matvecsubr( X, R, ipar )
R = R - B
call pcondlsubr( T1V, R, ipar )
residual = normfun( HUTI_NDIM, T1V, 1 ) / precrhsnorm
case (HUTI_XDIFF_NORM)
T1V = MATMUL( work(:,V_ind+1-1:V_ind+m-1), Y(1:m) )
residual = normfun( HUTI_NDIM, T1V, 1 )
case (HUTI_USUPPLIED_STOPC)
residual = stopcfun( X, B, R, ipar, dpar )
case default
call pcondrsubr( T1V, X, ipar )
call matvecsubr( T1V, R, ipar )
T1V = R - B
call pcondlsubr( R, T1V, ipar )
residual = normfun( HUTI_NDIM, R, 1 )
end select
work(m+1,S_ind) = normfun( HUTI_NDIM, R, 1 )
!
! Print debugging output if desired
!
if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
if ( mod(iter_count, HUTI_DBUGLVL) .eq. 0 ) then
WRITE (*, '(A, I8, E11.4)') ' gmresz:',iter_count, residual
end if
end if
if ( residual .lt. HUTI_TOLERANCE ) then
HUTI_INFO = HUTI_CONVERGENCE
go to 1000
end if
IF( residual /= residual .OR. residual > HUTI_MAXTOLERANCE ) THEN
HUTI_INFO = HUTI_DIVERGENCE
GOTO 1000
END IF
!
! Return back to the iteration loop (without initialization)
!
iter_count = iter_count + 1
if ( iter_count .gt. HUTI_MAXIT ) then
HUTI_INFO = HUTI_MAXITER
go to 1000
end if
go to 300
!
! This is where we exit last time (after enough iterations or breakdown)
!
1000 continue
if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
WRITE (*, '(A, I8, E11.4)') ' gmresz:', iter_count, residual
end if
HUTI_ITERS = iter_count
call pcondrsubr( T1V, X, ipar )
X = T1V
return
! End of execution
!*********************************************************************
end subroutine huti_cgmressolv
!*************************************************************************
!*************************************************************************
!*************************************************************************
! Double complex version
!*************************************************************************
!*************************************************************************
subroutine huti_zgmressolv ( ndim, wrkdim, xvec, rhsvec, &
ipar, dpar, work, matvecsubr, pcondlsubr, &
pcondrsubr, dotprodfun, normfun, stopcfun )
use huti_interfaces
implicit none
procedure( mv_iface_z ), pointer :: matvecsubr
procedure( pc_iface_z ), pointer :: pcondlsubr
procedure( pc_iface_z ), pointer :: pcondrsubr
procedure( dotp_iface_z ), pointer :: dotprodfun
procedure( norm_iface_z ), pointer :: normfun
procedure( stopc_iface_z ), pointer :: stopcfun
! Parameters
integer :: ndim, wrkdim
double complex, dimension(ndim) :: xvec, rhsvec
integer, dimension(HUTI_IPAR_DFLTSIZE) :: ipar
double precision, dimension(HUTI_DPAR_DFLTSIZE) :: dpar
double complex, dimension(ndim,wrkdim) :: work
! Local variables
integer :: iter_count
double precision :: residual, rhsnorm, precrhsnorm
double precision :: bnrm, alpha, beta
INTEGER :: reit, info, i, j, k, l, m, n
double complex :: temp, temp2, error, gamma
! Local arrays
double complex, DIMENSION(HUTI_GMRES_RESTART+1,HUTI_GMRES_RESTART+1) :: H
double complex, &
DIMENSION((HUTI_GMRES_RESTART+1)*(HUTI_GMRES_RESTART+1)) :: HLU
double complex, DIMENSION(HUTI_GMRES_RESTART+1) :: CS, SN, Y
!
! End of variable declarations
!*********************************************************************
!*********************************************************************
! The actual GMRES begins here (look the pseudo code in the
! "Templates..."-book, page 20)
!
! First the initialization part
!
iter_count = 1
bnrm = normfun( HUTI_NDIM, B, 1 )
! Norms of right-hand side vector are used in convergence tests
if ( HUTI_STOPC .eq. HUTI_TRESID_SCALED_BYB .or. &
HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYB ) then
rhsnorm = bnrm
end if
if ( HUTI_STOPC .eq. HUTI_PRESID_SCALED_BYPRECB ) then
call pcondlsubr( T1V, B, ipar )
precrhsnorm = normfun( HUTI_NDIM, T1V, 1 )
end if
! The following applies for all matrix operations in this solver
HUTI_EXTOP_MATTYPE = HUTI_MAT_NOTTRPSED
! Generate vector X if needed
if ( HUTI_INITIALX .eq. HUTI_RANDOMX ) then
call huti_zrandvec ( X, ipar )
else if ( HUTI_INITIALX .ne. HUTI_USERSUPPLIEDX ) then
X = 1
end if
call pcondrsubr( T1V, X, ipar )
call matvecsubr( T1V, R, ipar )
T1V = B - R
call pcondlsubr( R, T1V, ipar )
m = HUTI_GMRES_RESTART
work(:,V_ind+1-1:V_ind+m+1-1) = 0
H = 0
CS = 0
SN = 0
VTMP = 0
work(1,VTMP_ind) = 1.0
!
! This is where the loop starts (that is we continue from here after
! the first iteration)
!
300 continue
call pcondrsubr( T1V, X, ipar )
call matvecsubr( T1V, R, ipar )
T1V = B - R
call pcondlsubr( R, T1V, ipar )
alpha = normfun( HUTI_NDIM, R, 1 )
if ( alpha .eq. 0 ) then
HUTI_INFO = HUTI_GMRES_ALPHA
go to 1000
end if
work(:,V_ind+1-1) = R / alpha
S = alpha * VTMP
!
! Construct orthonormal
!
DO i = 1, m
call pcondrsubr( W, work(:,V_ind+i-1), ipar )
call matvecsubr( W, T1V, ipar )
call pcondlsubr( W, T1V, ipar )
DO k = 1, i
H(k,i) = dotprodfun( HUTI_NDIM, W, 1, work(:,V_ind+k-1), 1 )
W = W - H(k,i) * work(:,V_ind+k-1)
END DO
beta = normfun( HUTI_NDIM, W, 1 )
if ( beta .eq. 0 ) then
HUTI_INFO = HUTI_GMRES_BETA
go to 1000
end if
H(i+1,i) = beta
work(:,V_ind+i+1-1) = W / H(i+1, i)
!
! Compute the Givens rotation
!
DO k = 1, i-1
temp = CS(k) * H(k,i) + SN(k) * H(k+1,i)
H(k+1,i) = -1 * SN(k) * H(k,i) + CS(k) * H(k+1,i)
H(k,i) = temp
END DO
IF ( H(i+1,i) .eq. 0 ) THEN
CS(i) = 1; SN(i) = 0
ELSE
IF ( abs( H(i+1,i) ) .gt. abs( H(i,i) ) ) THEN
temp2 = H(i,i) / H(i+1,i)
SN(i) = 1 / sqrt( 1 + ( temp2 * temp2 ))
CS(i) = temp2 * SN(i)
ELSE
temp2 = H(i+1,i) / H(i,i)
CS(i) = 1 / sqrt( 1 + ( temp2 * temp2 ))
SN(i) = temp2 * CS(i)
END IF
END IF
temp = CS(i) * work(i,S_ind)
work(i+1,S_ind) = -1 * SN(i) * work(i,S_ind)
work(i,S_ind) = temp
H(i,i) = ( CS(i) * H(i,i) ) + ( SN(i) * H(i+1,i) )
H(i+1,i) = 0
error = abs( work(i+1,S_ind) ) / bnrm
IF ( REAL( error ) .lt. HUTI_TOLERANCE ) THEN
HLU = 0; j = 1
do k = 1, i
do l = 1, i
HLU(j) = H(l,k)
j = j + 1
end do
end do
call huti_zlusolve ( i, HLU, Y, work(:,S_ind) )
X = X + MATMUL( work(:,V_ind+1-1:V_ind+i-1), Y(1:i) )
EXIT
END IF
END DO
IF ( REAL( error ) .lt. HUTI_TOLERANCE ) THEN
GOTO 500
END IF
HLU = 0; j = 1
do k = 1, m
do l = 1, m
HLU(j) = H(l,k)
j = j + 1
end do
end do
call huti_zlusolve ( m, HLU, Y, work(:,S_ind) )
X = X + MATMUL( work(:,V_ind+1-1:V_ind+m-1), Y(1:m) )
500 CONTINUE
!
! Check the convergence
!
select case (HUTI_STOPC)
case (HUTI_TRUERESIDUAL)
call pcondrsubr( T1V, X, ipar )
call matvecsubr( T1V, R, ipar )
T1V = B - R
call pcondlsubr( R, T1V, ipar )
residual = normfun( HUTI_NDIM, R, 1 )
case (HUTI_TRESID_SCALED_BYB)
call pcondrsubr( T1V, X, ipar )
call matvecsubr( T1V, R, ipar )
T1V = B - R
call pcondlsubr( R, T1V, ipar )
residual = normfun( HUTI_NDIM, R, 1 ) / rhsnorm
case (HUTI_PSEUDORESIDUAL)
call matvecsubr( X, R, ipar )
R = R - B
call pcondlsubr( T1V, R, ipar )
residual = normfun( HUTI_NDIM, T1V, 1 )
case (HUTI_PRESID_SCALED_BYB)
call matvecsubr( X, R, ipar )
R = R - B
call pcondlsubr( T1V, R, ipar )
residual = normfun( HUTI_NDIM, T1V, 1 ) / rhsnorm
case (HUTI_PRESID_SCALED_BYPRECB)
call matvecsubr( X, R, ipar )
R = R - B
call pcondlsubr( T1V, R, ipar )
residual = normfun( HUTI_NDIM, T1V, 1 ) / precrhsnorm
case (HUTI_XDIFF_NORM)
T1V = MATMUL( work(:,V_ind+1-1:V_ind+m-1), Y(1:m) )
residual = normfun( HUTI_NDIM, T1V, 1 )
case (HUTI_USUPPLIED_STOPC)
residual = stopcfun( X, B, R, ipar, dpar )
case default
call pcondrsubr( T1V, X, ipar )
call matvecsubr( T1V, R, ipar )
T1V = R - B
call pcondlsubr( R, T1V, ipar )
residual = normfun( HUTI_NDIM, R, 1 )
end select
work(m+1,S_ind) = normfun( HUTI_NDIM, R, 1 )
!
! Print debugging output if desired
!
if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
if ( mod(iter_count, HUTI_DBUGLVL) .eq. 0 ) then
WRITE (*, '(A, I8, E11.4)') ' gmresz:',iter_count, residual
end if
end if
if ( residual .lt. HUTI_TOLERANCE ) then
HUTI_INFO = HUTI_CONVERGENCE
go to 1000
end if
IF( residual /= residual .OR. residual > HUTI_MAXTOLERANCE ) THEN
HUTI_INFO = HUTI_DIVERGENCE
GOTO 1000
END IF
!
! Return back to the iteration loop (without initialization)
!
iter_count = iter_count + 1
if ( iter_count .gt. HUTI_MAXIT ) then
HUTI_INFO = HUTI_MAXITER
go to 1000
end if
go to 300
!
! This is where we exit last time (after enough iterations or breakdown)
!
1000 continue
if ( HUTI_DBUGLVL .ne. HUTI_NO_DEBUG ) then
WRITE (*, '(A,I8, E11.4)') ' gmresz:',iter_count, residual
end if
HUTI_ITERS = iter_count
call pcondrsubr( T1V, X, ipar )
X = T1V
return
! End of execution
!*********************************************************************
end subroutine huti_zgmressolv
!*************************************************************************
end module huti_gmres
Computing file changes ...
