fdurastante
/
psblas3
mirror of https://github.com/sfilippone/psblas3.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Krylov iterations:

Browse Source 
Centralized check convergence routines. Still incomplete for RGMRES.
psblas3-type-indexed
Salvatore Filippone 17 years ago
parent 274375c40b
commit c304258426
7 changed files with 680 additions and 660 deletions
Show Stats Download Patch File Download Diff File

114
krylov/psb_dbicg.f90
Unescape Escape View File

@ -115,16 +115,15 @@ subroutine psb_dbicg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
real(kind(1.d0)), pointer :: ww(:), q(:),&
& r(:), p(:), zt(:), pt(:), z(:), rt(:),qt(:)
integer :: int_err(5)
integer ::litmax, naux, mglob, it, itrace_,&
integer ::itmax_, naux, mglob, it, itrace_,&
& np,me, n_row, n_col, istop_, err_act
integer :: debug_level, debug_unit
logical, parameter :: exchange=.true., noexchange=.false.
integer, parameter :: irmax = 8
integer :: itx, isvch, ictxt
real(kind(1.d0)) :: alpha, beta, rho, rho_old, rni, xni, bni, ani,&
& sigma,bn2
real(kind(1.d0)) :: errnum, errden
character(len=20) :: name,ch_err
real(kind(1.d0)) :: alpha, beta, rho, rho_old, sigma
type(psb_itconv_type) :: stopdat
character(len=20) :: name,ch_err
character(len=*), parameter :: methdname='BiCG'
info = 0
@ -202,9 +201,9 @@ subroutine psb_dbicg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
ww => wwrk(:,9)
if (present(itmax)) then
litmax = itmax
itmax_ = itmax
else
litmax = 1000
itmax_ = 1000
endif
if (present(itrace)) then
@ -215,27 +214,18 @@ subroutine psb_dbicg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
itx = 0
if (istop_ == 1) then
ani = psb_spnrmi(a,desc_a,info)
bni = psb_geamax(b,desc_a,info)
else if (istop_ == 2) then
bn2 = psb_genrm2(b,desc_a,info)
endif
errnum = dzero
errden = done
if(info /= 0) then
info=4011
err=info
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call psb_init_conv(istop_,itrace_,a,b,eps,desc_a,stopdat,info)
if (info /= 0) Then
call psb_errpush(4011,name)
goto 9999
End If
restart: do
!!$
!!$ r0 = b-ax0
!!$
if (itx >= litmax) exit restart
if (itx >= itmax_) exit restart
it = 0
call psb_geaxpby(done,b,dzero,r,desc_a,info)
if (info == 0) call psb_spmm(-done,a,x,done,r,desc_a,info,work=aux)
@ -249,48 +239,23 @@ subroutine psb_dbicg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
end if
rho = dzero
if (debug_level >= psb_debug_ext_)&
& write(debug_unit,*) me,' ',trim(name),'on entry to amax: b: ',size(b)
if (istop_ == 1) then
rni = psb_geamax(r,desc_a,info)
xni = psb_geamax(x,desc_a,info)
else if (istop_ == 2) then
rni = psb_genrm2(r,desc_a,info)
endif
if(info /= 0) then
info=4011
call psb_errpush(info,name)
goto 9999
end if
if (istop_ == 1) then
xni = psb_geamax(x,desc_a,info)
errnum = rni
errden = (ani*xni+bni)
else if (istop_ == 2) then
errnum = rni
errden = bn2
endif
if(info /= 0) then
info=4011
call psb_errpush(info,name)
! Perhaps we already satisfy the convergence criterion...
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (info /= 0) Then
call psb_errpush(4011,name)
goto 9999
end if
if (errnum <= eps*errden) exit restart
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,eps)
End If
iteration: do
it = it + 1
itx = itx + 1
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),'iteration: ',itx
call psb_precaply(prec,r,z,desc_a,info,work=aux)
call psb_precaply(prec,rt,zt,desc_a,info,trans='t',work=aux)
if (info == 0) call psb_precaply(prec,rt,zt,desc_a,info,trans='t',work=aux)
rho_old = rho
rho = psb_gedot(rt,z,desc_a,info)
@ -330,43 +295,16 @@ subroutine psb_dbicg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
call psb_geaxpby(-alpha,q,done,r,desc_a,info)
call psb_geaxpby(-alpha,qt,done,rt,desc_a,info)
if (istop_ == 1) then
rni = psb_geamax(r,desc_a,info)
xni = psb_geamax(x,desc_a,info)
errnum = rni
errden = (ani*xni+bni)
else if (istop_ == 2) then
rni = psb_genrm2(r,desc_a,info)
errnum = rni
errden = bn2
endif
if (errnum <= eps*errden) exit restart
if (itx >= litmax) exit restart
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,eps)
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (info /= 0) Then
call psb_errpush(4011,name)
goto 9999
End If
end do iteration
end do restart
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,1,errnum,errden,eps)
if (present(err)) then
if (errden /= dzero) then
err = errnum/errden
else
err = errnum
end if
end if
if (present(iter)) iter = itx
if (errnum > eps*errden) &
& call end_log(methdname,me,itx,errnum,errden,eps)
call psb_end_conv(methdname,itx,desc_a,stopdat,info,err,iter)
deallocate(aux, stat=info)
if (info == 0) call psb_gefree(wwrk,desc_a,info)

107
krylov/psb_dcg.f90
Unescape Escape View File

@ -114,18 +114,19 @@ subroutine psb_dcg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
!!$ Local data
real(kind(1.d0)), allocatable, target :: aux(:), wwrk(:,:)
real(kind(1.d0)), pointer :: q(:), p(:), r(:), z(:), w(:)
real(kind(1.d0)) ::alpha, beta, rho, rho_old, rni, xni, bni, ani,bn2,&
& sigma
integer :: litmax, istop_, naux, mglob, it, itx, itrace_,&
real(kind(1.d0)) :: alpha, beta, rho, rho_old, sigma
integer :: itmax_, istop_, naux, mglob, it, itx, itrace_,&
& np,me, n_col, isvch, ictxt, n_row,err_act, int_err(5)
logical, parameter :: exchange=.true., noexchange=.false.
real(kind(1.d0)) :: errnum, errden
integer :: debug_level, debug_unit
type(psb_itconv_type) :: stopdat
character(len=20) :: name
character(len=*), parameter :: methdname='CG'
info = 0
name = 'psb_dcg'
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = psb_cd_get_context(desc_a)
@ -141,30 +142,12 @@ subroutine psb_dcg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
else
istop_ = 1
endif
!
! ISTOP_ = 1: Normwise backward error, infinity norm
! ISTOP_ = 2: ||r||/||b|| norm 2
!
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
write(0,*) 'psb_cg: invalid istop',istop_
info=5001
int_err(1)=istop_
err=info
call psb_errpush(info,name,i_err=int_err)
goto 9999
endif
call psb_chkvect(mglob,1,size(x,1),1,1,desc_a,info)
if(info /= 0) then
info=4010
call psb_errpush(info,name,a_err='psb_chkvect on X')
goto 9999
end if
call psb_chkvect(mglob,1,size(b,1),1,1,desc_a,info)
if (info == 0) call psb_chkvect(mglob,1,size(b,1),1,1,desc_a,info)
if(info /= 0) then
info=4010
call psb_errpush(info,name,a_err='psb_chkvect on B')
call psb_errpush(info,name,a_err='psb_chkvect on X/B')
goto 9999
end if
@ -186,9 +169,9 @@ subroutine psb_dcg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
if (present(itmax)) then
litmax = itmax
itmax_ = itmax
else
litmax = 1000
itmax_ = 1000
endif
if (present(itrace)) then
@ -206,10 +189,11 @@ subroutine psb_dcg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
!!$
!!$ r0 = b-Ax0
!!$
if (itx>= litmax) exit restart
if (itx>= itmax_) exit restart
it = 0
call psb_geaxpby(done,b,dzero,r,desc_a,info)
call psb_spmm(-done,a,x,done,r,desc_a,info,work=aux)
if (info == 0) call psb_spmm(-done,a,x,done,r,desc_a,info,work=aux)
if (info /= 0) then
info=4011
call psb_errpush(info,name)
@ -217,22 +201,15 @@ subroutine psb_dcg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
end if
rho = dzero
if (istop_ == 1) then
ani = psb_spnrmi(a,desc_a,info)
bni = psb_geamax(b,desc_a,info)
else if (istop_ == 2) then
bn2 = psb_genrm2(b,desc_a,info)
endif
errnum = dzero
errden = done
if (info /= 0) then
info=4011
call psb_errpush(info,name)
call psb_init_conv(istop_,itrace_,a,b,eps,desc_a,stopdat,info)
if (info /= 0) Then
call psb_errpush(4011,name)
goto 9999
end if
End If
iteration: do
it = it + 1
itx = itx + 1
@ -244,7 +221,9 @@ subroutine psb_dcg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
call psb_geaxpby(done,z,dzero,p,desc_a,info)
else
if (rho_old==dzero) then
write(0,*) 'CG Iteration breakdown'
if (debug_level >= psb_debug_ext_)&
& write(debug_unit,*) me,' ',trim(name),&
& ': CG Iteration breakdown rho'
exit iteration
endif
beta = rho/rho_old
@ -254,7 +233,9 @@ subroutine psb_dcg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
call psb_spmm(done,a,p,dzero,q,desc_a,info,work=aux)
sigma = psb_gedot(p,q,desc_a,info)
if (sigma==dzero) then
write(0,*) 'CG Iteration breakdown'
if (debug_level >= psb_debug_ext_)&
& write(debug_unit,*) me,' ',trim(name),&
& ': CG Iteration breakdown sigma'
exit iteration
endif
@ -262,42 +243,16 @@ subroutine psb_dcg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
call psb_geaxpby(alpha,p,done,x,desc_a,info)
call psb_geaxpby(-alpha,q,done,r,desc_a,info)
if (istop_ == 1) then
rni = psb_geamax(r,desc_a,info)
xni = psb_geamax(x,desc_a,info)
errnum = rni
errden = (ani*xni+bni)
else if (istop_ == 2) then
rni = psb_genrm2(r,desc_a,info)
errnum = rni
errden = bn2
endif
if (errnum <= eps*errden) exit restart
if (itx>= litmax) exit restart
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,eps)
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (info /= 0) Then
call psb_errpush(4011,name)
goto 9999
End If
end do iteration
end do restart
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,1,errnum,errden,eps)
if (present(err)) then
if (errden /= dzero) then
err = errnum/errden
else
err = errnum
end if
end if
if (present(iter)) iter = itx
if (errnum > eps*errden) &
& call end_log(methdname,me,itx,errnum,errden,eps)
call psb_end_conv(methdname,itx,desc_a,stopdat,info,err,iter)
call psb_gefree(wwrk,desc_a,info)
if (info /= 0) then

179
krylov/psb_dcgs.f90
Unescape Escape View File

@ -114,15 +114,12 @@ Subroutine psb_dcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
Real(Kind(1.d0)), allocatable, target :: aux(:),wwrk(:,:)
Real(Kind(1.d0)), Pointer :: ww(:), q(:),&
& r(:), p(:), v(:), s(:), z(:), f(:), rt(:),qt(:),uv(:)
Integer :: litmax, naux, mglob, it, itrace_,int_err(5),&
Integer :: itmax_, naux, mglob, it, itrace_,int_err(5),&
& np,me, n_row, n_col,istop_, err_act
Logical, Parameter :: exchange=.True., noexchange=.False.
Integer, Parameter :: irmax = 8
Integer :: itx, isvch, ictxt
integer :: debug_level, debug_unit
Real(Kind(1.d0)) :: alpha, beta, rho, rho_old, rni, xni, bni, ani,bn2,&
& sigma
real(kind(1.d0)) :: errnum, errden
Real(Kind(1.d0)) :: alpha, beta, rho, rho_old, sigma
type(psb_itconv_type) :: stopdat
character(len=20) :: name
character(len=*), parameter :: methdname='CGS'
@ -146,29 +143,12 @@ Subroutine psb_dcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
Else
istop_ = 1
Endif
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
info=5001
int_err=istop_
err=info
call psb_errpush(info,name,i_err=int_err)
goto 9999
endif
call psb_chkvect(mglob,1,size(x,1),1,1,desc_a,info)
if(info /= 0) then
info=4010
call psb_errpush(info,name,a_err='psb_chkvect on X')
goto 9999
end if
call psb_chkvect(mglob,1,size(b,1),1,1,desc_a,info)
if (info == 0) call psb_chkvect(mglob,1,size(b,1),1,1,desc_a,info)
if(info /= 0) then
info=4010
call psb_errpush(info,name,a_err='psb_chkvect on B')
call psb_errpush(info,name,a_err='psb_chkvect on X/B')
goto 9999
end if
@ -196,9 +176,9 @@ Subroutine psb_dcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
If (Present(itmax)) Then
litmax = itmax
itmax_ = itmax
Else
litmax = 1000
itmax_ = 1000
Endif
If (Present(itrace)) Then
@ -212,69 +192,47 @@ Subroutine psb_dcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
itx = 0
if (istop_ == 1) then
ani = psb_spnrmi(a,desc_a,info)
bni = psb_geamax(b,desc_a,info)
else if (istop_ == 2) then
bn2 = psb_genrm2(b,desc_a,info)
endif
errnum = dzero
errden = done
if(info/=0)then
info=4011
call psb_errpush(info,name)
call psb_init_conv(istop_,itrace_,a,b,eps,desc_a,stopdat,info)
if (info /= 0) Then
call psb_errpush(4011,name)
goto 9999
end if
End If
restart: Do
!!$
!!$ r0 = b-ax0
!!$
if (itx >= litmax) exit restart
if (itx >= itmax_) exit restart
it = 0
call psb_geaxpby(done,b,dzero,r,desc_a,info)
call psb_spmm(-done,a,x,done,r,desc_a,info,work=aux)
call psb_geaxpby(done,r,dzero,rt,desc_a,info)
if(info/=0)then
if (info == 0) call psb_spmm(-done,a,x,done,r,desc_a,info,work=aux)
if (info == 0) call psb_geaxpby(done,r,dzero,rt,desc_a,info)
if (info/=0) then
info=4011
call psb_errpush(info,name)
goto 9999
end if
! Perhaps we already satisfy the convergence criterion...
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (info /= 0) Then
call psb_errpush(4011,name)
goto 9999
End If
rho = dzero
If (debug_level >= psb_debug_ext_)&
& write(debug_unit,*) me,' ',trim(name),&
& ' on entry to amax: b: ',Size(b)
if (istop_ == 1) then
rni = psb_geamax(r,desc_a,info)
xni = psb_geamax(x,desc_a,info)
errnum = rni
errden = (ani*xni+bni)
else if (istop_ == 2) then
rni = psb_genrm2(r,desc_a,info)
errnum = rni
errden = bn2
endif
if(info/=0)then
info=4011
call psb_errpush(info,name)
goto 9999
end if
if (errnum <= eps*errden) then
exit restart
end if
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,eps)
iteration: do
it = it + 1
itx = itx + 1
If (debug_level >= psb_debug_ext_) &
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),'iteration: ',itx
rho_old = rho
rho = psb_gedot(rt,r,desc_a,info)
if (rho==dzero) then
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
@ -284,20 +242,25 @@ Subroutine psb_dcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
if (it==1) then
call psb_geaxpby(done,r,dzero,uv,desc_a,info)
call psb_geaxpby(done,r,dzero,p,desc_a,info)
if (info == 0) call psb_geaxpby(done,r,dzero,p,desc_a,info)
else
beta = (rho/rho_old)
call psb_geaxpby(done,r,dzero,uv,desc_a,info)
call psb_geaxpby(beta,q,done,uv,desc_a,info)
call psb_geaxpby(done,q,beta,p,desc_a,info)
call psb_geaxpby(done,uv,beta,p,desc_a,info)
if (info == 0) call psb_geaxpby(beta,q,done,uv,desc_a,info)
if (info == 0) call psb_geaxpby(done,q,beta,p,desc_a,info)
if (info == 0) call psb_geaxpby(done,uv,beta,p,desc_a,info)
end if
call psb_precaply(prec,p,f,desc_a,info,work=aux)
if (info == 0) call psb_precaply(prec,p,f,desc_a,info,work=aux)
call psb_spmm(done,a,f,dzero,v,desc_a,info,&
if (info == 0) call psb_spmm(done,a,f,dzero,v,desc_a,info,&
& work=aux)
if (info /= 0) then
call psb_errpush(4010,name,a_err='First loop part ')
goto 9999
end if
sigma = psb_gedot(rt,v,desc_a,info)
if (sigma==dzero) then
if (debug_level >= psb_debug_ext_) &
@ -308,61 +271,41 @@ Subroutine psb_dcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
alpha = rho/sigma
call psb_geaxpby(done,uv,dzero,q,desc_a,info)
call psb_geaxpby(-alpha,v,done,q,desc_a,info)
call psb_geaxpby(done,uv,dzero,s,desc_a,info)
call psb_geaxpby(done,q,done,s,desc_a,info)
if (info == 0) call psb_geaxpby(done,uv,dzero,q,desc_a,info)
if (info == 0) call psb_geaxpby(-alpha,v,done,q,desc_a,info)
if (info == 0) call psb_geaxpby(done,uv,dzero,s,desc_a,info)
if (info == 0) call psb_geaxpby(done,q,done,s,desc_a,info)
call psb_precaply(prec,s,z,desc_a,info,work=aux)
if (info == 0) call psb_precaply(prec,s,z,desc_a,info,work=aux)
call psb_geaxpby(alpha,z,done,x,desc_a,info)
if (info == 0) call psb_geaxpby(alpha,z,done,x,desc_a,info)
call psb_spmm(done,a,z,dzero,qt,desc_a,info,&
if (info == 0) call psb_spmm(done,a,z,dzero,qt,desc_a,info,&
& work=aux)
call psb_geaxpby(-alpha,qt,done,r,desc_a,info)
if (info == 0) call psb_geaxpby(-alpha,qt,done,r,desc_a,info)
if (istop_ == 1) then
rni = psb_geamax(r,desc_a,info)
xni = psb_geamax(x,desc_a,info)
errnum = rni
errden = (ani*xni+bni)
else if (istop_ == 2) then
rni = psb_genrm2(r,desc_a,info)
errnum = rni
errden = bn2
endif
if (errnum <= eps*errden) exit restart
if (itx >= litmax) exit restart
if (info /= 0) then
call psb_errpush(4010,name,a_err='X update ')
goto 9999
end if
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,eps)
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (info /= 0) Then
call psb_errpush(4011,name)
goto 9999
End If
end do iteration
end do restart
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,1,errnum,errden,eps)
if (present(iter)) iter = itx
if (present(err)) then
if (errden /= dzero) then
err = errnum/errden
else
err = errnum
end if
end if
if (errnum > eps*errden) &
& call end_log(methdname,me,itx,errnum,errden,eps)
call psb_end_conv(methdname,itx,desc_a,stopdat,info,err,iter)
deallocate(aux,stat=info)
if (info == 0) call psb_gefree(wwrk,desc_a,info)
if(info/=0) then
call psb_errpush(info,name)
goto 9999
if (info /= 0) then
call psb_errpush(info,name)
goto 9999
end if
! restore external global coherence behaviour
@ -373,8 +316,8 @@ Subroutine psb_dcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
call psb_error()
return
end if
return

172
krylov/psb_dcgstab.F90
Unescape Escape View File

@ -113,16 +113,16 @@ Subroutine psb_dcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
Real(Kind(1.d0)), allocatable, target :: aux(:),wwrk(:,:)
Real(Kind(1.d0)), Pointer :: q(:),&
& r(:), p(:), v(:), s(:), t(:), z(:), f(:)
Integer :: litmax, naux, mglob, it,itrace_,&
Integer :: itmax_, naux, mglob, it,itrace_,&
& np,me, n_row, n_col
integer :: debug_level, debug_unit
Logical, Parameter :: exchange=.True., noexchange=.False., debug1 = .False.
Integer, Parameter :: irmax = 8
Integer :: itx, isvch, ictxt, err_act, int_err(5)
Integer :: itx, isvch, ictxt, err_act
Integer :: istop_
Real(Kind(1.d0)) :: alpha, beta, rho, rho_old, rni, xni, bni, ani,&
& sigma, omega, tau, rn0, bn2
real(kind(1.d0)) :: errnum, errden
Real(Kind(1.d0)) :: alpha, beta, rho, rho_old, sigma, omega, tau
type(psb_itconv_type) :: stopdat
#ifdef MPE_KRYLOV
Integer istpb, istpe, ifctb, ifcte, imerr, irank, icomm,immb,imme
Integer mpe_log_get_event_number,mpe_Describe_state,mpe_log_event
@ -174,13 +174,6 @@ Subroutine psb_dcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
imerr = MPE_Log_event( istpb, 0, "st CGSTAB" )
#endif
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
info=5001
int_err(1)=istop_
err=info
call psb_errpush(info,name,i_err=int_err)
goto 9999
endif
call psb_chkvect(mglob,1,size(x,1),1,1,desc_a,info)
if(info /= 0) then
@ -215,9 +208,9 @@ Subroutine psb_dcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
Z => WWRK(:,8)
If (Present(itmax)) Then
litmax = itmax
itmax_ = itmax
Else
litmax = 1000
itmax_ = 1000
Endif
If (Present(itrace)) Then
@ -230,99 +223,51 @@ Subroutine psb_dcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
call psb_set_coher(ictxt,isvch)
itx = 0
If (istop_ == 1) Then
ani = psb_spnrmi(a,desc_a,info)
bni = psb_geamax(b,desc_a,info)
Else If (istop_ == 2) Then
bn2 = psb_genrm2(b,desc_a,info)
Endif
errnum = dzero
errden = done
call psb_init_conv(istop_,itrace_,a,b,eps,desc_a,stopdat,info)
if (info /= 0) Then
info=4011
call psb_errpush(info,name)
call psb_errpush(4011,name)
goto 9999
End If
restart: Do
!!$
!!$ r0 = b-Ax0
!!$
if (itx >= litmax) exit restart
if (itx >= itmax_) exit restart
it = 0
call psb_geaxpby(done,b,dzero,r,desc_a,info)
#ifdef MPE_KRYLOV
imerr = MPE_Log_event( immb, 0, "st SPMM" )
#endif
call psb_spmm(-done,a,x,done,r,desc_a,info,work=aux)
if (info == 0) call psb_spmm(-done,a,x,done,r,desc_a,info,work=aux)
#ifdef MPE_KRYLOV
imerr = MPE_Log_event( imme, 0, "ed SPMM" )
#endif
call psb_geaxpby(done,r,dzero,q,desc_a,info)
if (info /= 0) then
info=4011
call psb_errpush(info,name)
goto 9999
end if
rho = dzero
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' On entry to AMAX: B: ',Size(b)
!
! Must always provide norm of R into RNI below for first check on
! residual
!
if (istop_ == 1) then
rni = psb_geamax(r,desc_a,info)
xni = psb_geamax(x,desc_a,info)
else if (istop_ == 2) then
rni = psb_genrm2(r,desc_a,info)
endif
if (info == 0) call psb_geaxpby(done,r,dzero,q,desc_a,info)
if (info /= 0) then
info=4011
call psb_errpush(info,name)
info=4010
call psb_errpush(info,name,a_err='Init residual')
goto 9999
end if
if (itx == 0) then
rn0 = rni
end if
if (rn0 == dzero ) then
errnum = dzero
errden = done
exit restart
end if
! Perhaps we already satisfy the convergence criterion...
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (info /= 0) Then
call psb_errpush(4011,name)
goto 9999
End If
if (istop_ == 1) then
xni = psb_geamax(x,desc_a,info)
errnum = rni
errden = (ani*xni+bni)
else if (istop_ == 2) then
errnum = rni
errden = bn2
endif
if (info /= 0) then
info=4011
call psb_errpush(info,name)
goto 9999
end if
if (errnum <= eps*errden) exit restart
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,1,errnum,errden,eps)
rho = dzero
iteration: Do
it = it + 1
itx = itx + 1
if (debug_level >= psb_debug_ext_)&
& write(debug_unit,*) me,' ',trim(name),&
& ' Iteration: ',itx
rho_old = rho
rho = psb_gedot(q,r,desc_a,info)
rho = psb_gedot(q,r,desc_a,info)
if (rho==dzero) then
if (debug_level >= psb_debug_ext_) &
@ -363,11 +308,8 @@ Subroutine psb_dcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
alpha = rho/sigma
call psb_geaxpby(done,r,dzero,s,desc_a,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_geaxpby')
goto 9999
end if
call psb_geaxpby(-alpha,v,done,s,desc_a,info)
if (info == 0) call psb_geaxpby(-alpha,v,done,s,desc_a,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_geaxpby')
goto 9999
@ -377,25 +319,19 @@ Subroutine psb_dcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
imerr = MPE_Log_event( ifctb, 0, "st PREC" )
#endif
call psb_precaply(prec,s,z,desc_a,info,work=aux)
#ifdef MPE_KRYLOV
imerr = MPE_Log_event( ifcte, 0, "ed PREC" )
#endif
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_precaply')
goto 9999
end if
#ifdef MPE_KRYLOV
imerr = MPE_Log_event( ifcte, 0, "ed PREC" )
imerr = MPE_Log_event( immb, 0, "st SPMM" )
#endif
Call psb_spmm(done,a,z,dzero,t,desc_a,info,&
if (info == 0) Call psb_spmm(done,a,z,dzero,t,desc_a,info,&
& work=aux)
#ifdef MPE_KRYLOV
imerr = MPE_Log_event( imme, 0, "ed SPMM" )
#endif
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_spmm')
call psb_errpush(4010,name,a_err='precaply/spmm')
goto 9999
end if
@ -407,7 +343,7 @@ Subroutine psb_dcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
exit iteration
endif
tau = psb_gedot(t,s,desc_a,info)
tau = psb_gedot(t,s,desc_a,info)
omega = tau/sigma
if (omega==dzero) then
@ -418,42 +354,24 @@ Subroutine psb_dcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
endif
call psb_geaxpby(alpha,f,done,x,desc_a,info)
call psb_geaxpby(omega,z,done,x,desc_a,info)
call psb_geaxpby(done,s,dzero,r,desc_a,info)
call psb_geaxpby(-omega,t,done,r,desc_a,info)
if (info == 0) call psb_geaxpby(omega,z,done,x,desc_a,info)
if (info == 0) call psb_geaxpby(done,s,dzero,r,desc_a,info)
if (info == 0) call psb_geaxpby(-omega,t,done,r,desc_a,info)
if (info /= 0) Then
call psb_errpush(4010,name,a_err='X/R update ')
goto 9999
End If
if (istop_ == 1) then
rni = psb_geamax(r,desc_a,info)
xni = psb_geamax(x,desc_a,info)
errnum = rni
errden = (ani*xni+bni)
else if (istop_ == 2) then
rni = psb_genrm2(r,desc_a,info)
errnum = rni
errden = bn2
endif
if (errnum <= eps*errden) exit restart
if (itx >= litmax) exit restart
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,eps)
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (info /= 0) Then
call psb_errpush(4011,name)
goto 9999
End If
end do iteration
end do restart
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,1,errnum,errden,eps)
if (present(err)) then
if (errden /= dzero) then
err = errnum/errden
else
err = errnum
end if
end if
if (present(iter)) iter = itx
if (errnum > eps*errden) &
& call end_log(methdname,me,itx,errnum,errden,eps)
call psb_end_conv(methdname,itx,desc_a,stopdat,info,err,iter)
deallocate(aux,stat=info)
if (info == 0) call psb_gefree(wwrk,desc_a,info)

403
krylov/psb_krylov_mod.f90
Unescape Escape View File

@ -34,6 +34,7 @@
!
Module psb_krylov_mod
public
interface psb_krylov
module procedure psb_dkrylov, psb_zkrylov
@ -192,6 +193,30 @@ Module psb_krylov_mod
end subroutine psb_zcgs
end interface
interface psb_init_conv
module procedure psb_d_init_conv, psb_z_init_conv
end interface
interface psb_check_conv
module procedure psb_d_check_conv, psb_z_check_conv
end interface
interface psb_end_conv
module procedure psb_end_conv
end interface
integer, parameter, private :: bni_=1, rni_=2, ani_=3, xni_=4, bn2_=5, xn2_=6
integer, parameter, private :: errnum_=7, errden_=8, eps_=9, rn2_=10
integer, parameter, private :: stopc_=1, trace_=2
integer, parameter, private :: ivsz_=16
type psb_itconv_type
private
integer :: controls(ivsz_)
real(kind(1.d0)) :: values(ivsz_)
end type psb_itconv_type
contains
! Subroutine: psb_dkrylov
!
@ -424,17 +449,24 @@ contains
end subroutine psb_zkrylov
subroutine log_conv(methdname,me,itx,itrace,errnum,errden,eps)
use psb_base_mod
character(len=*), intent(in) :: methdname
integer, intent(in) :: me, itx, itrace
real(kind(1.d0)), intent(in) :: errnum, errden, eps
character(len=*), parameter :: fmt='(a,i4,3(2x,es10.4))'
if ((mod(itx,itrace)==0).and.(me == 0)) then
write(*,fmt) methdname//': ',itx,errnum,eps*errden
if (me == 0) then
if (errden > dzero ) then
write(*,fmt) methdname//': ',itx,errnum/errden,eps
else
write(*,fmt) methdname//': ',itx,errnum,eps
end if
Endif
end subroutine log_conv
subroutine end_log(methdname,me,itx,errnum,errden,eps)
use psb_base_mod
character(len=*), intent(in) :: methdname
integer, intent(in) :: me, itx
real(kind(1.d0)), intent(in) :: errnum, errden, eps
@ -444,12 +476,375 @@ contains
if (me==0) then
write(*,fmt) methdname//' failed to converge to ',eps,&
& ' in ',ITX,' iterations. '
write(*,fmt1) 'Last iteration convergence indicators: ',&
& errnum,eps*errden,errnum/errden
if (errden > dzero) then
write(*,fmt1) 'Last iteration convergence indicator: ',&
& errnum/errden
else
write(*,fmt1) 'Last iteration convergence indicator: ',&
& errnum/errden
endif
end if
end subroutine end_log
subroutine psb_d_init_conv(stopc,trace,a,b,eps,desc_a,stopdat,info)
use psb_base_mod
implicit none
integer, intent(in) :: stopc, trace
type(psb_dspmat_type), intent(in) :: a
real(kind(1.d0)), intent(in) :: b(:), eps
type(psb_desc_type), intent(in) :: desc_a
type(psb_itconv_type) :: stopdat
integer, intent(out) :: info
integer :: ictxt, me, np, err_act
character(len=20) :: name
info = 0
name = 'psb_init_conv'
call psb_erractionsave(err_act)
ictxt=psb_cd_get_context(desc_a)
call psb_info(ictxt, me, np)
stopdat%controls(:) = 0
stopdat%values(:) = 0.0d0
stopdat%controls(stopc_) = stopc
stopdat%controls(trace_) = trace
select case(stopdat%controls(stopc_))
case (1)
stopdat%values(ani_) = psb_spnrmi(a,desc_a,info)
if (info == 0) stopdat%values(bni_) = psb_geamax(b,desc_a,info)
case (2)
stopdat%values(bn2_) = psb_genrm2(b,desc_a,info)
case default
info=5001
call psb_errpush(info,name,i_err=(/stopc,0,0,0,0/))
goto 9999
end select
if (info /= 0) then
call psb_errpush(4001,name,a_err="Init conv check data")
goto 9999
end if
stopdat%values(eps_) = eps
stopdat%values(errnum_) = dzero
stopdat%values(errden_) = done
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error(ictxt)
return
end if
end subroutine psb_d_init_conv
subroutine psb_z_init_conv(stopc,trace,a,b,eps,desc_a,stopdat,info)
use psb_base_mod
implicit none
integer, intent(in) :: stopc, trace
type(psb_zspmat_type), intent(in) :: a
complex(kind(1.d0)), intent(in) :: b(:)
real(kind(1.d0)), intent(in) :: eps
type(psb_desc_type), intent(in) :: desc_a
type(psb_itconv_type) :: stopdat
integer, intent(out) :: info
integer :: ictxt, me, np, err_act
character(len=20) :: name
info = 0
name = 'psb_init_conv'
call psb_erractionsave(err_act)
ictxt=psb_cd_get_context(desc_a)
call psb_info(ictxt, me, np)
stopdat%controls(:) = 0
stopdat%values(:) = 0.0d0
stopdat%controls(stopc_) = stopc
stopdat%controls(trace_) = trace
select case(stopdat%controls(stopc_))
case (1)
stopdat%values(ani_) = psb_spnrmi(a,desc_a,info)
if (info == 0) stopdat%values(bni_) = psb_geamax(b,desc_a,info)
case (2)
stopdat%values(bn2_) = psb_genrm2(b,desc_a,info)
case default
info=5001
call psb_errpush(info,name,i_err=(/stopc,0,0,0,0/))
goto 9999
end select
if (info /= 0) then
call psb_errpush(4001,name,a_err="Init conv check data")
goto 9999
end if
stopdat%values(eps_) = eps
stopdat%values(errnum_) = dzero
stopdat%values(errden_) = done
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error(ictxt)
return
end if
end subroutine psb_z_init_conv
function psb_d_check_conv(methdname,it,x,r,desc_a,stopdat,info)
use psb_base_mod
implicit none
character(len=*), intent(in) :: methdname
integer, intent(in) :: it
real(kind(1.d0)), intent(in) :: x(:), r(:)
type(psb_desc_type), intent(in) :: desc_a
type(psb_itconv_type) :: stopdat
logical :: psb_d_check_conv
integer, intent(out) :: info
integer :: ictxt, me, np, err_act
real(kind(1.d0)) :: errnum, errden
character(len=*), parameter :: fmt='(a,i4,3(2x,es10.4))'
character(len=20) :: name
info = 0
name = 'psb_check_conv'
call psb_erractionsave(err_act)
ictxt = psb_cd_get_context(desc_a)
call psb_info(ictxt,me,np)
psb_d_check_conv = .false.
select case(stopdat%controls(stopc_))
case(1)
stopdat%values(rni_) = psb_geamax(r,desc_a,info)
if (info == 0) stopdat%values(xni_) = psb_geamax(x,desc_a,info)
stopdat%values(errnum_) = stopdat%values(rni_)
stopdat%values(errden_) = &
& (stopdat%values(ani_)*stopdat%values(xni_)+stopdat%values(bni_))
case(2)
stopdat%values(rn2_) = psb_genrm2(r,desc_a,info)
stopdat%values(errnum_) = stopdat%values(rn2_)
stopdat%values(errden_) = stopdat%values(bn2_)
case default
info=4001
call psb_errpush(info,name,a_err="Control data in stopdat messed up!")
goto 9999
end select
if (info /= 0) then
info=4011
call psb_errpush(info,name)
goto 9999
end if
if (stopdat%values(errden_) == dzero) then
psb_d_check_conv = (stopdat%values(errnum_) <= stopdat%values(eps_))
else
psb_d_check_conv = &
& (stopdat%values(errnum_) <= stopdat%values(eps_)*stopdat%values(errden_))
end if
if (((stopdat%controls(trace_) > 0).and.(mod(it,stopdat%controls(trace_))==0))&
& .or.psb_d_check_conv) then
if (me == 0) then
if (stopdat%values(errden_) > dzero ) then
write(*,fmt) trim(methdname)//': ',it,&
& stopdat%values(errnum_)/stopdat%values(errden_),stopdat%values(eps_)
else
write(*,fmt) trim(methdname)//': ',it,&
& stopdat%values(errnum_),stopdat%values(eps_)
end if
Endif
end if
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error(ictxt)
return
end if
end function psb_d_check_conv
function psb_z_check_conv(methdname,it,x,r,desc_a,stopdat,info)
use psb_base_mod
implicit none
character(len=*), intent(in) :: methdname
integer, intent(in) :: it
complex(kind(1.d0)), intent(in) :: x(:), r(:)
type(psb_desc_type), intent(in) :: desc_a
type(psb_itconv_type) :: stopdat
logical :: psb_z_check_conv
integer, intent(out) :: info
integer :: ictxt, me, np, err_act
real(kind(1.d0)) :: errnum, errden
character(len=*), parameter :: fmt='(a,i4,3(2x,es10.4))'
character(len=20) :: name
info = 0
name = 'psb_check_conv'
call psb_erractionsave(err_act)
ictxt = psb_cd_get_context(desc_a)
call psb_info(ictxt,me,np)
psb_z_check_conv = .false.
select case(stopdat%controls(stopc_))
case(1)
stopdat%values(rni_) = psb_geamax(r,desc_a,info)
if (info == 0) stopdat%values(xni_) = psb_geamax(x,desc_a,info)
stopdat%values(errnum_) = stopdat%values(rni_)
stopdat%values(errden_) = &
& (stopdat%values(ani_)*stopdat%values(xni_)+stopdat%values(bni_))
case(2)
stopdat%values(rn2_) = psb_genrm2(r,desc_a,info)
stopdat%values(errnum_) = stopdat%values(rn2_)
stopdat%values(errden_) = stopdat%values(bn2_)
case default
info=4001
call psb_errpush(info,name,a_err="Control data in stopdat messed up!")
goto 9999
end select
if (info /= 0) then
info=4011
call psb_errpush(info,name)
goto 9999
end if
if (stopdat%values(errden_) == dzero) then
psb_z_check_conv = (stopdat%values(errnum_) <= stopdat%values(eps_))
else
psb_z_check_conv = &
& (stopdat%values(errnum_) <= stopdat%values(eps_)*stopdat%values(errden_))
end if
if (((stopdat%controls(trace_) > 0).and.(mod(it,stopdat%controls(trace_))==0))&
& .or.psb_z_check_conv) then
if (me == 0) then
if (stopdat%values(errden_) > dzero ) then
write(*,fmt) trim(methdname)//': ',it,&
& stopdat%values(errnum_)/stopdat%values(errden_),stopdat%values(eps_)
else
write(*,fmt) trim(methdname)//': ',it,&
& stopdat%values(errnum_),stopdat%values(eps_)
end if
Endif
end if
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error(ictxt)
return
end if
end function psb_z_check_conv
subroutine psb_end_conv(methdname,it,desc_a,stopdat,info,err,iter)
use psb_base_mod
implicit none
character(len=*), intent(in) :: methdname
integer, intent(in) :: it
type(psb_desc_type), intent(in) :: desc_a
type(psb_itconv_type) :: stopdat
integer, intent(out) :: info
real(kind(1.d0)), optional, intent(out) :: err
integer, optional, intent(out) :: iter
integer :: ictxt, me, np, err_act
real(kind(1.d0)) :: errnum, errden, eps
character(len=*), parameter :: fmt='(a,2x,es10.4,1x,a,1x,i4,1x,a)'
character(len=*), parameter :: fmt1='(a,3(2x,es10.4))'
character(len=20) :: name
info = 0
name = 'psb_end_conv'
call psb_erractionsave(err_act)
ictxt = psb_cd_get_context(desc_a)
call psb_info(ictxt,me,np)
if (present(iter)) iter = it
errnum = stopdat%values(errnum_)
errden = stopdat%values(errden_)
eps = stopdat%values(eps_)
if (errden == dzero) then
if (errnum > eps) then
if (me==0) then
write(*,fmt) methdname//' failed to converge to ',eps,&
& ' in ',it,' iterations. '
write(*,fmt1) 'Last iteration convergence indicator: ',&
& errnum
end if
end if
if (present(err)) err=errnum
else
if (errnum > eps) then
if (me==0) then
write(*,fmt) methdname//' failed to converge to ',eps,&
& ' in ',it,' iterations. '
write(*,fmt1) 'Last iteration convergence indicator: ',&
& errnum/errden
end if
endif
if (present(err)) err=errnum/errden
end if
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error(ictxt)
return
end if
end subroutine psb_end_conv
end module psb_krylov_mod

194
krylov/psb_zcgs.f90
Unescape Escape View File

@ -112,16 +112,13 @@ Subroutine psb_zcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
Complex(Kind(1.d0)), allocatable, target :: aux(:),wwrk(:,:)
Complex(Kind(1.d0)), Pointer :: ww(:), q(:),&
& r(:), p(:), v(:), s(:), z(:), f(:), rt(:),qt(:),uv(:)
Integer :: litmax, naux, mglob, it, itrace_,int_err(5),&
Integer :: itmax_, naux, mglob, it, itrace_,int_err(5),&
& np,me, n_row, n_col,istop_, err_act
Logical, Parameter :: exchange=.True., noexchange=.False.
Integer, Parameter :: irmax = 8
Integer :: itx, isvch, ictxt
integer :: debug_level, debug_unit
Real(Kind(1.d0)) :: rni, xni, bni, ani,bn2
complex(Kind(1.d0)) :: alpha, beta, rho, rho_old, sigma
real(kind(1.d0)) :: errnum, errden
character(len=20) :: name
complex(Kind(1.d0)) :: alpha, beta, rho, rho_old, sigma
type(psb_itconv_type) :: stopdat
character(len=20) :: name
character(len=*), parameter :: methdname='CGS'
info = 0
@ -144,29 +141,12 @@ Subroutine psb_zcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
Else
istop_ = 1
Endif
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
info=5001
int_err=istop_
err=info
call psb_errpush(info,name,i_err=int_err)
goto 9999
endif
call psb_chkvect(mglob,1,size(x,1),1,1,desc_a,info)
if(info /= 0) then
info=4010
call psb_errpush(info,name,a_err='psb_chkvect on X')
goto 9999
end if
call psb_chkvect(mglob,1,size(b,1),1,1,desc_a,info)
if (info == 0) call psb_chkvect(mglob,1,size(b,1),1,1,desc_a,info)
if(info /= 0) then
info=4010
call psb_errpush(info,name,a_err='psb_chkvect on B')
call psb_errpush(info,name,a_err='psb_chkvect on X/B')
goto 9999
end if
@ -194,83 +174,63 @@ Subroutine psb_zcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
If (Present(itmax)) Then
litmax = itmax
itmax_ = itmax
Else
litmax = 1000
itmax_ = 1000
Endif
if (present(itrace)) then
If (Present(itrace)) Then
itrace_ = itrace
else
Else
itrace_ = 0
end if
End If
! Ensure global coherence for convergence checks.
call psb_set_coher(ictxt,isvch)
itx = 0
if (istop_ == 1) then
ani = psb_spnrmi(a,desc_a,info)
bni = psb_geamax(b,desc_a,info)
else if (istop_ == 2) then
bn2 = psb_genrm2(b,desc_a,info)
endif
errnum = dzero
errden = done
if(info/=0)then
info=4011
call psb_errpush(info,name)
goto 9999
end if
call psb_init_conv(istop_,itrace_,a,b,eps,desc_a,stopdat,info)
if (info /= 0) Then
call psb_errpush(4011,name)
goto 9999
End If
restart: Do
!!$
!!$ r0 = b-ax0
!!$
If (itx >= litmax) Exit restart
if (itx >= itmax_) exit restart
it = 0
Call psb_geaxpby(zone,b,zzero,r,desc_a,info)
Call psb_spmm(-zone,a,x,zone,r,desc_a,info,work=aux)
Call psb_geaxpby(zone,r,zzero,rt,desc_a,info)
if(info/=0)then
info=4011
call psb_errpush(info,name)
goto 9999
call psb_geaxpby(zone,b,zzero,r,desc_a,info)
if (info == 0) call psb_spmm(-zone,a,x,zone,r,desc_a,info,work=aux)
if (info == 0) call psb_geaxpby(zone,r,zzero,rt,desc_a,info)
if (info/=0) then
info=4011
call psb_errpush(info,name)
goto 9999
end if
rho = zzero
If (debug_level >= psb_debug_ext_)&
& write(debug_unit,*) me,' ',trim(name),&
& ' on entry to amax: b: ',Size(b)
if (istop_ == 1) then
rni = psb_geamax(r,desc_a,info)
xni = psb_geamax(x,desc_a,info)
errnum = rni
errden = (ani*xni+bni)
else if (istop_ == 2) then
rni = psb_genrm2(r,desc_a,info)
errnum = rni
errden = bn2
endif
if(info/=0)then
info=4011
call psb_errpush(info,name)
! Perhaps we already satisfy the convergence criterion...
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (info /= 0) Then
call psb_errpush(4011,name)
goto 9999
end if
End If
if (errnum <= eps*errden) exit restart
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,eps)
rho = zzero
iteration: Do
iteration: do
it = it + 1
itx = itx + 1
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),'iteration: ',itx
rho_old = rho
rho = psb_gedot(rt,r,desc_a,info)
if (rho==zzero) then
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
@ -280,20 +240,25 @@ Subroutine psb_zcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
if (it==1) then
call psb_geaxpby(zone,r,zzero,uv,desc_a,info)
call psb_geaxpby(zone,r,zzero,p,desc_a,info)
if (info == 0) call psb_geaxpby(zone,r,zzero,p,desc_a,info)
else
beta = (rho/rho_old)
call psb_geaxpby(zone,r,zzero,uv,desc_a,info)
call psb_geaxpby(beta,q,zone,uv,desc_a,info)
call psb_geaxpby(zone,q,beta,p,desc_a,info)
call psb_geaxpby(zone,uv,beta,p,desc_a,info)
if (info == 0) call psb_geaxpby(beta,q,zone,uv,desc_a,info)
if (info == 0) call psb_geaxpby(zone,q,beta,p,desc_a,info)
if (info == 0) call psb_geaxpby(zone,uv,beta,p,desc_a,info)
end if
call psb_precaply(prec,p,f,desc_a,info,work=aux)
if (info == 0) call psb_precaply(prec,p,f,desc_a,info,work=aux)
call psb_spmm(zone,a,f,zzero,v,desc_a,info,&
if (info == 0) call psb_spmm(zone,a,f,zzero,v,desc_a,info,&
& work=aux)
if (info /= 0) then
call psb_errpush(4010,name,a_err='First loop part ')
goto 9999
end if
sigma = psb_gedot(rt,v,desc_a,info)
if (sigma==zzero) then
if (debug_level >= psb_debug_ext_) &
@ -301,56 +266,38 @@ Subroutine psb_zcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
& ' iteration breakdown s1', sigma
exit iteration
endif
alpha = rho/sigma
call psb_geaxpby(zone,uv,zzero,q,desc_a,info)
call psb_geaxpby(-alpha,v,zone,q,desc_a,info)
call psb_geaxpby(zone,uv,zzero,s,desc_a,info)
call psb_geaxpby(zone,q,zone,s,desc_a,info)
call psb_precaply(prec,s,z,desc_a,info,work=aux)
if (info == 0) call psb_geaxpby(zone,uv,zzero,q,desc_a,info)
if (info == 0) call psb_geaxpby(-alpha,v,zone,q,desc_a,info)
if (info == 0) call psb_geaxpby(zone,uv,zzero,s,desc_a,info)
if (info == 0) call psb_geaxpby(zone,q,zone,s,desc_a,info)
if (info == 0) call psb_precaply(prec,s,z,desc_a,info,work=aux)
call psb_geaxpby(alpha,z,zone,x,desc_a,info)
if (info == 0) call psb_geaxpby(alpha,z,zone,x,desc_a,info)
call psb_spmm(zone,a,z,zzero,qt,desc_a,info,&
if (info == 0) call psb_spmm(zone,a,z,zzero,qt,desc_a,info,&
& work=aux)
if (info == 0) call psb_geaxpby(-alpha,qt,zone,r,desc_a,info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='X update ')
goto 9999
end if
call psb_geaxpby(-alpha,qt,zone,r,desc_a,info)
if (istop_ == 1) then
rni = psb_geamax(r,desc_a,info)
xni = psb_geamax(x,desc_a,info)
errnum = rni
errden = (ani*xni+bni)
else if (istop_ == 2) then
rni = psb_genrm2(r,desc_a,info)
errnum = rni
errden = bn2
endif
if (errnum <= eps*errden) exit restart
if (itx >= litmax) exit restart
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,eps)
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (info /= 0) Then
call psb_errpush(4011,name)
goto 9999
End If
end do iteration
end do restart
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,1,errnum,errden,eps)
if (present(err)) then
if (errden /= dzero) then
err = errnum/errden
else
err = errnum
end if
end if
if (present(iter)) iter = itx
if (errnum > eps*errden) &
& call end_log(methdname,me,itx,errnum,errden,eps)
call psb_end_conv(methdname,itx,desc_a,stopdat,info,err,iter)
deallocate(aux,stat=info)
if (info == 0) call psb_gefree(wwrk,desc_a,info)
@ -361,7 +308,6 @@ Subroutine psb_zcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
! restore external global coherence behaviour
call psb_restore_coher(ictxt,isvch)
call psb_erractionrestore(err_act)
return

171
krylov/psb_zcgstab.f90
Unescape Escape View File

@ -115,13 +115,10 @@ subroutine psb_zcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
Integer :: litmax, naux, mglob, it,itrace_,&
& np,me, n_row, n_col
integer :: debug_level, debug_unit
Logical, Parameter :: exchange=.True., noexchange=.False., debug1 = .False.
Integer, Parameter :: irmax = 8
Integer :: itx, isvch, ictxt, err_act, int_err(5)
Integer :: itx, isvch, ictxt, err_act
Integer :: istop_
complex(Kind(1.d0)) :: alpha, beta, rho, rho_old, sigma, omega, tau
Real(Kind(1.d0)) :: rni, xni, bni, ani, rn0, bn2
real(kind(1.d0)) :: errnum, errden
type(psb_itconv_type) :: stopdat
!!$ Integer istpb, istpe, ifctb, ifcte, imerr, irank, icomm,immb,imme
!!$ Integer mpe_log_get_event_number,mpe_Describe_state,mpe_log_event
character(len=20) :: name
@ -146,18 +143,6 @@ subroutine psb_zcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
Else
istop_ = 1
Endif
!
! ISTOP = 1: Normwise backward error, infinity norm
! ISTOP = 2: ||r||/||b|| norm 2
!
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
info=5001
int_err(1)=istop_
err=info
call psb_errpush(info,name,i_err=int_err)
goto 9999
endif
call psb_chkvect(mglob,1,size(x,1),1,1,desc_a,info)
if(info /= 0) then
@ -208,18 +193,14 @@ subroutine psb_zcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
itx = 0
If (istop_ == 1) Then
ani = psb_spnrmi(a,desc_a,info)
bni = psb_geamax(b,desc_a,info)
Else If (istop_ == 2) Then
bn2 = psb_genrm2(b,desc_a,info)
Endif
call psb_init_conv(istop_,itrace_,a,b,eps,desc_a,stopdat,info)
if (info /= 0) Then
info=4011
call psb_errpush(info,name)
call psb_errpush(4011,name)
goto 9999
End If
restart: Do
!!$
!!$ r0 = b-Ax0
@ -227,8 +208,8 @@ subroutine psb_zcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
if (itx >= litmax) exit restart
it = 0
call psb_geaxpby(zone,b,zzero,r,desc_a,info)
call psb_spmm(-zone,a,x,zone,r,desc_a,info,work=aux)
call psb_geaxpby(zone,r,zzero,q,desc_a,info)
if (info == 0) call psb_spmm(-zone,a,x,zone,r,desc_a,info,work=aux)
if (info == 0) call psb_geaxpby(zone,r,zzero,q,desc_a,info)
if (info /= 0) then
info=4011
call psb_errpush(info,name)
@ -240,46 +221,13 @@ subroutine psb_zcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
& write(debug_unit,*) me,' ',trim(name),&
& ' On entry to AMAX: B: ',Size(b)
!
! Must always provide norm of R into RNI below for first check on
! residual
!
if (istop_ == 1) then
rni = psb_geamax(r,desc_a,info)
xni = psb_geamax(x,desc_a,info)
else if (istop_ == 2) then
rni = psb_genrm2(r,desc_a,info)
endif
errnum = dzero
errden = done
if (info /= 0) then
info=4011
call psb_errpush(info,name)
goto 9999
end if
if (itx == 0) then
rn0 = rni
end if
if (rn0 == 0.d0 ) exit restart
if (istop_ == 1) then
xni = psb_geamax(x,desc_a,info)
errnum = rni
errden = (ani*xni+bni)
else if (istop_ == 2) then
errnum = rni
errden = bn2
endif
if (info /= 0) then
info=4011
call psb_errpush(info,name)
goto 9999
end if
if (errnum <= eps*errden) exit restart
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,eps)
! Perhaps we already satisfy the convergence criterion...
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (info /= 0) Then
call psb_errpush(4011,name)
goto 9999
End If
iteration: Do
it = it + 1
@ -287,8 +235,9 @@ subroutine psb_zcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
If (debug_level >= psb_debug_ext_)&
& write(debug_unit,*) me,' ',trim(name),&
& ' Iteration: ',itx
rho_old = rho
rho = psb_gedot(q,r,desc_a,info)
rho = psb_gedot(q,r,desc_a,info)
if (rho==zzero) then
if (debug_level >= psb_debug_ext_) &
@ -302,47 +251,40 @@ subroutine psb_zcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
else
beta = (rho/rho_old)*(alpha/omega)
call psb_geaxpby(-omega,v,zone,p,desc_a,info)
call psb_geaxpby(zone,r,beta,p,desc_a,info)
if (info == 0) call psb_geaxpby(zone,r,beta,p,desc_a,info)
end if
call psb_precaply(prec,p,f,desc_a,info,work=aux)
if (info == 0) call psb_precaply(prec,p,f,desc_a,info,work=aux)
call psb_spmm(zone,a,f,zzero,v,desc_a,info,&
if (info == 0) call psb_spmm(zone,a,f,zzero,v,desc_a,info,&
& work=aux)
sigma = psb_gedot(q,v,desc_a,info)
if (info == 0) sigma = psb_gedot(q,v,desc_a,info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='First step')
goto 9999
end if
if (sigma==zzero) then
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Iteration breakdown S1', sigma
exit iteration
endif
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' SIGMA:',sigma
alpha = rho/sigma
call psb_geaxpby(zone,r,zzero,s,desc_a,info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_geaxpby')
goto 9999
end if
call psb_geaxpby(-alpha,v,zone,s,desc_a,info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_geaxpby')
goto 9999
end if
call psb_precaply(prec,s,z,desc_a,info,work=aux)
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_precaply')
goto 9999
end if
call psb_spmm(zone,a,z,zzero,t,desc_a,info,&
call psb_geaxpby(zone,r,zzero,s,desc_a,info)
if (info == 0) call psb_geaxpby(-alpha,v,zone,s,desc_a,info)
if (info == 0) call psb_precaply(prec,s,z,desc_a,info,work=aux)
if (info == 0) call psb_spmm(zone,a,z,zzero,t,desc_a,info,&
& work=aux)
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_spmm')
call psb_errpush(4010,name,a_err='Second step ')
goto 9999
end if
@ -354,7 +296,7 @@ subroutine psb_zcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
exit iteration
endif
tau = psb_gedot(t,s,desc_a,info)
tau = psb_gedot(t,s,desc_a,info)
omega = tau/sigma
if (omega==zzero) then
@ -364,43 +306,26 @@ subroutine psb_zcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
exit iteration
endif
call psb_geaxpby(alpha,f,zone,x,desc_a,info)
call psb_geaxpby(omega,z,zone,x,desc_a,info)
call psb_geaxpby(zone,s,zzero,r,desc_a,info)
call psb_geaxpby(-omega,t,zone,r,desc_a,info)
if (istop_ == 1) then
rni = psb_geamax(r,desc_a,info)
xni = psb_geamax(x,desc_a,info)
errnum = rni
errden = (ani*xni+bni)
else if (istop_ == 2) then
rni = psb_genrm2(r,desc_a,info)
errnum = rni
errden = bn2
endif
if (errnum <= eps*errden) exit restart
if (itx >= litmax) exit restart
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,eps)
if (info == 0) call psb_geaxpby(alpha,f,zone,x,desc_a,info)
if (info == 0) call psb_geaxpby(omega,z,zone,x,desc_a,info)
if (info == 0) call psb_geaxpby(zone,s,zzero,r,desc_a,info)
if (info == 0) call psb_geaxpby(-omega,t,zone,r,desc_a,info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='X update ')
goto 9999
end if
if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
if (info /= 0) Then
call psb_errpush(4011,name)
goto 9999
End If
end do iteration
end do restart
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,1,errnum,errden,eps)
if (present(err)) then
if (errden /= dzero) then
err = errnum/errden
else
err = errnum
end if
end if
if (present(iter)) iter = itx
if (errnum > eps*errden) &
& call end_log(methdname,me,itx,errnum,errden,eps)
call psb_end_conv(methdname,itx,desc_a,stopdat,info,err,iter)
deallocate(aux,stat=info)
if (info == 0) call psb_gefree(wwrk,desc_a,info)

Write Preview
Loading…
Cancel
Save