fdurastante/psblas3
1
0
Fork 0
mirror of https://github.com/sfilippone/psblas3.git synced 2026-10-18 03:42:09 +00:00
Code Issues Packages Projects Releases Wiki Activity

Cleaned source code.

Browse Source
This commit is contained in:
Salvatore Filippone
2008-01-11 20:48:47 +00:00
parent bc980ed5d6
commit 7f758de51c
10 changed files with 609 additions and 710 deletions
Download Patch File Download Diff File Expand all files Collapse all files
krylov/psb_dbicg.f90
+28 -36
View File
@@ -96,6 +96,7 @@
subroutine psb_dbicg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
use psb_base_mod
use psb_prec_mod
use psb_krylov_mod, psb_protect_name => psb_dbicg
implicit none
!!$ parameters
@@ -124,6 +125,7 @@ subroutine psb_dbicg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
& sigma,bn2
real(kind(1.d0)) :: errnum, errden
character(len=20) :: name,ch_err
character(len=*), parameter :: methdname='BiCG'
info = 0
name = 'psb_dbicg'
@@ -181,7 +183,7 @@ subroutine psb_dbicg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
allocate(aux(naux),stat=info)
if (info == 0) call psb_geall(wwrk,desc_a,info,n=9)
if (info == 0) call psb_geasb(wwrk,desc_a,info)
if(info.ne.0) then
if(info /= 0) then
info=4011
ch_err='psb_asb'
err=info
@@ -222,7 +224,7 @@ subroutine psb_dbicg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
errnum = dzero
errden = done
if(info.ne.0) then
if(info /= 0) then
info=4011
err=info
call psb_errpush(info,name,a_err=ch_err)
@@ -233,14 +235,14 @@ subroutine psb_dbicg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
!!$
!!$ r0 = b-ax0
!!$
if (itx.ge.litmax) exit restart
if (itx >= litmax) 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)
if (debug_level >= psb_debug_ext_)&
& write(debug_unit,*) me,' ',trim(name),' Done spmm',info
if (info == 0) call psb_geaxpby(done,r,dzero,rt,desc_a,info)
if(info.ne.0) then
if(info /= 0) then
info=4011
call psb_errpush(info,name)
goto 9999
@@ -255,7 +257,7 @@ subroutine psb_dbicg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
else if (istop_ == 2) then
rni = psb_genrm2(r,desc_a,info)
endif
if(info.ne.0) then
if(info /= 0) then
info=4011
call psb_errpush(info,name)
goto 9999
@@ -270,20 +272,16 @@ subroutine psb_dbicg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
errden = bn2
endif
if(info.ne.0) then
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) then
if ((mod(itx,itrace_)==0).and.(me == 0))&
& write(*,'(a,i4,3(2x,es10.4))') 'bicg: ',itx,errnum,eps*errden
end If
if (errnum <= eps*errden) exit restart
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,eps)
iteration: do
it = it + 1
@@ -343,21 +341,19 @@ subroutine psb_dbicg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
errnum = rni
errden = bn2
endif
If (errnum <= eps*errden) Then
exit restart
end if
if (itx.ge.litmax) exit restart
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 (itrace_ > 0) then
if ((mod(itx,itrace_)==0).and.(me == 0))&
& write(*,'(a,i4,3(2x,es10.4))') 'bicg: ',itx,errnum,eps*errden
end If
end do iteration
end do restart
If (itrace_ > 0) then
if (me == 0) write(*,'(a,i4,3(2x,es10.4))') 'bicg: ',itx,errnum,eps*errden
end If
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,1,errnum,errden,eps)
if (present(err)) then
if (errden /= dzero) then
@@ -368,29 +364,25 @@ subroutine psb_dbicg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
end if
if (present(iter)) iter = itx
If ((errnum > eps*errden).and.(me==0)) Then
write(debug_unit,*) 'bicg failed to converge to ',eps*errden,&
& ' in ',itx,' iterations '
end if
if (errnum > eps*errden) &
& call end_log(methdname,me,itx,errnum,errden,eps)
deallocate(aux)
call psb_gefree(wwrk,desc_a,info)
! restore external global coherence behaviour
call psb_restore_coher(ictxt,isvch)
if(info/=0) then
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
end if
! restore external global coherence behaviour
call psb_restore_coher(ictxt,isvch)
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
krylov/psb_dcg.f90
+26 -31
View File
@@ -94,9 +94,10 @@
! estimate of) residual.
!
!
Subroutine psb_dcg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
subroutine psb_dcg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
use psb_base_mod
use psb_prec_mod
use psb_krylov_mod, psb_protect_name => psb_dcg
implicit none
!!$ Parameters
@@ -119,7 +120,8 @@ Subroutine psb_dcg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
& 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
character(len=20) :: name
character(len=20) :: name
character(len=*), parameter :: methdname='CG'
info = 0
name = 'psb_dcg'
@@ -170,7 +172,7 @@ Subroutine psb_dcg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
allocate(aux(naux), stat=info)
if (info == 0) call psb_geall(wwrk,desc_a,info,n=5)
if (info == 0) call psb_geasb(wwrk,desc_a,info)
if (info.ne.0) then
if (info /= 0) then
info=4011
call psb_errpush(info,name)
goto 9999
@@ -208,7 +210,7 @@ Subroutine psb_dcg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
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.ne.0) then
if (info /= 0) then
info=4011
call psb_errpush(info,name)
goto 9999
@@ -223,7 +225,7 @@ Subroutine psb_dcg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
endif
errnum = dzero
errden = done
if (info.ne.0) then
if (info /= 0) then
info=4011
call psb_errpush(info,name)
goto 9999
@@ -234,7 +236,7 @@ Subroutine psb_dcg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
it = it + 1
itx = itx + 1
Call psb_precaply(prec,r,z,desc_a,info,work=aux)
call psb_precaply(prec,r,z,desc_a,info,work=aux)
rho_old = rho
rho = psb_gedot(r,z,desc_a,info)
@@ -271,54 +273,47 @@ Subroutine psb_dcg(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
errnum = rni
errden = bn2
endif
if (errnum <= eps*errden) then
exit restart
end if
if (errnum <= eps*errden) exit restart
if (itx>= litmax) exit restart
If (itrace_ > 0) then
if ((mod(itx,itrace_)==0).and.(me == 0))&
& write(*,'(a,i4,3(2x,es10.4))') 'cg: ',itx,errnum,eps*errden
end If
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,eps)
end do iteration
end do restart
If (itrace_ > 0) then
if (me == 0) write(*,'(a,i4,3(2x,es10.4))') 'cg: ',itx,errnum,eps*errden
end If
If (Present(err)) then
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).and.(me==0)) Then
write(0,*) 'CG Failed to converge to ',eps*errden,&
& ' in ',litmax,' iterations '
info=itx
end if
deallocate(aux)
if (present(iter)) iter = itx
if (errnum > eps*errden) &
& call end_log(methdname,me,itx,errnum,errden,eps)
call psb_gefree(wwrk,desc_a,info)
! restore external global coherence behaviour
call psb_restore_coher(ictxt,isvch)
if (info.ne.0) then
if (info /= 0) then
call psb_errpush(info,name)
goto 9999
end if
! restore external global coherence behaviour
call psb_restore_coher(ictxt,isvch)
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
krylov/psb_dcgs.f90
+62 -73
View File
@@ -96,6 +96,7 @@
Subroutine psb_dcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
use psb_base_mod
use psb_prec_mod
use psb_krylov_mod, psb_protect_name => psb_dcgs
implicit none
!!$ parameters
@@ -119,10 +120,11 @@ Subroutine psb_dcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
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,&
Real(Kind(1.d0)) :: alpha, beta, rho, rho_old, rni, xni, bni, ani,bn2,&
& sigma
real(kind(1.d0)) :: errnum, errden
character(len=20) :: name
real(kind(1.d0)) :: errnum, errden
character(len=20) :: name
character(len=*), parameter :: methdname='CGS'
info = 0
name = 'psb_dcgs'
@@ -174,7 +176,7 @@ Subroutine psb_dcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
Allocate(aux(naux),stat=info)
if (info == 0) Call psb_geall(wwrk,desc_a,info,n=11)
if (info == 0) Call psb_geasb(wwrk,desc_a,info)
if (info.ne.0) Then
if (info /= 0) Then
info=4011
call psb_errpush(info,name)
goto 9999
@@ -228,11 +230,11 @@ Subroutine psb_dcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
!!$
!!$ r0 = b-ax0
!!$
If (itx.Ge.litmax) Exit restart
if (itx >= litmax) 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)
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
info=4011
call psb_errpush(info,name)
@@ -259,68 +261,66 @@ Subroutine psb_dcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
call psb_errpush(info,name)
goto 9999
end if
If (errnum <= eps*errden) Then
Exit restart
End If
If (itrace_ > 0) then
if ((mod(itx,itrace_)==0).and.(me == 0))&
& write(*,'(a,i4,3(2x,es10.4))') 'cgs: ',itx,errnum,eps*errden
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
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==dzero) Then
If (debug_level >= psb_debug_ext_) &
if (rho==dzero) then
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' iteration breakdown r',rho
Exit iteration
Endif
exit iteration
endif
If (it==1) Then
Call psb_geaxpby(done,r,dzero,uv,desc_a,info)
Call psb_geaxpby(done,r,dzero,p,desc_a,info)
Else
if (it==1) then
call psb_geaxpby(done,r,dzero,uv,desc_a,info)
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)
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)
end if
End If
call psb_precaply(prec,p,f,desc_a,info,work=aux)
Call psb_precaply(prec,p,f,desc_a,info,work=aux)
Call psb_spmm(done,a,f,dzero,v,desc_a,info,&
call psb_spmm(done,a,f,dzero,v,desc_a,info,&
& work=aux)
sigma = psb_gedot(rt,v,desc_a,info)
If (sigma==dzero) Then
If (debug_level >= psb_debug_ext_) &
if (sigma==dzero) then
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' iteration breakdown s1', sigma
Exit iteration
Endif
exit iteration
endif
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)
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)
Call psb_precaply(prec,s,z,desc_a,info,work=aux)
call psb_precaply(prec,s,z,desc_a,info,work=aux)
Call psb_geaxpby(alpha,z,done,x,desc_a,info)
call psb_geaxpby(alpha,z,done,x,desc_a,info)
Call psb_spmm(done,a,z,dzero,qt,desc_a,info,&
call psb_spmm(done,a,z,dzero,qt,desc_a,info,&
& work=aux)
Call psb_geaxpby(-alpha,qt,done,r,desc_a,info)
call psb_geaxpby(-alpha,qt,done,r,desc_a,info)
if (istop_ == 1) then
@@ -334,56 +334,45 @@ Subroutine psb_dcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
errden = bn2
endif
If (errnum <= eps*errden) Then
Exit restart
End If
If (itx.Ge.litmax) Exit restart
if (errnum <= eps*errden) exit restart
if (itx >= litmax) exit restart
If (itrace_ > 0) then
if ((mod(itx,itrace_)==0).and.(me == 0))&
& write(*,'(a,i4,3(2x,es10.4))') 'cgs: ',itx,errnum,eps*errden
end If
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,eps)
End Do iteration
End Do restart
end do iteration
end do restart
If (itrace_ > 0) then
if ((mod(itx,itrace_)==0).and.(me == 0))&
& write(*,'(a,i4,3(2x,es10.4))') 'cgs: ',itx,errnum,eps*errden
end If
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,1,errnum,errden,eps)
If (Present(iter)) iter = itx
If (Present(err)) then
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).and.(me==0)) Then
write(debug_unit,*) 'cgs failed to converge to ',eps*errden,&
& ' in ',itx,' iterations '
End If
Deallocate(aux)
Call psb_gefree(wwrk,desc_a,info)
! restore external global coherence behaviour
call psb_restore_coher(ictxt,isvch)
end if
if (errnum > eps*errden) &
& call end_log(methdname,me,itx,errnum,errden,eps)
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
end if
! restore external global coherence behaviour
call psb_restore_coher(ictxt,isvch)
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
krylov/psb_dcgstab.F90
+88 -112
View File
@@ -96,6 +96,7 @@
Subroutine psb_dcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
use psb_base_mod
use psb_prec_mod
use psb_krylov_mod, psb_protect_name => psb_dcgstab
Implicit None
!!$ parameters
Type(psb_dspmat_type), Intent(in) :: a
@@ -126,7 +127,8 @@ Subroutine psb_dcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
Integer istpb, istpe, ifctb, ifcte, imerr, irank, icomm,immb,imme
Integer mpe_log_get_event_number,mpe_Describe_state,mpe_log_event
#endif
character(len=20) :: name
character(len=20) :: name
character(len=*), parameter :: methdname='BiCGStab'
info = 0
name = 'psb_dcgstab'
@@ -247,25 +249,25 @@ Subroutine psb_dcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
!!$
!!$ r0 = b-Ax0
!!$
If (itx >= litmax) Exit restart
if (itx >= litmax) exit restart
it = 0
Call psb_geaxpby(done,b,dzero,r,desc_a,info)
call psb_geaxpby(done,b,dzero,r,desc_a,info)
#ifdef MPE_KRYLOV
imerr = MPE_Log_event( immb, 0, "st SPMM" )
imerr = MPE_Log_event( immb, 0, "st SPMM" )
#endif
Call psb_spmm(-done,a,x,done,r,desc_a,info,work=aux)
call psb_spmm(-done,a,x,done,r,desc_a,info,work=aux)
#ifdef MPE_KRYLOV
imerr = MPE_Log_event( imme, 0, "ed SPMM" )
imerr = MPE_Log_event( imme, 0, "ed SPMM" )
#endif
Call psb_geaxpby(done,r,dzero,q,desc_a,info)
if (info /= 0) Then
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
end if
rho = dzero
If (debug_level >= psb_debug_ext_) &
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' On entry to AMAX: B: ',Size(b)
@@ -273,113 +275,99 @@ Subroutine psb_dcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
! Must always provide norm of R into RNI below for first check on
! residual
!
If (istop_ == 1) Then
if (istop_ == 1) then
rni = psb_geamax(r,desc_a,info)
xni = psb_geamax(x,desc_a,info)
Else If (istop_ == 2) Then
else if (istop_ == 2) then
rni = psb_genrm2(r,desc_a,info)
Endif
if (info /= 0) Then
endif
if (info /= 0) then
info=4011
call psb_errpush(info,name)
goto 9999
End If
end if
If (itx == 0) Then
if (itx == 0) then
rn0 = rni
End If
If (rn0 == dzero ) Then
If (itrace_ > 0 ) Then
If (me == 0) Write(*,*) 'BiCGSTAB: ',itx,rn0
Endif
end if
if (rn0 == dzero ) then
errnum = dzero
errden = done
Exit restart
End If
exit restart
end if
If (istop_ == 1) Then
if (istop_ == 1) then
xni = psb_geamax(x,desc_a,info)
errnum = rni
errden = (ani*xni+bni)
Else If (istop_ == 2) Then
else if (istop_ == 2) then
errnum = rni
errden = bn2
Endif
if (info /= 0) Then
endif
if (info /= 0) then
info=4011
call psb_errpush(info,name)
goto 9999
End If
end if
If (errnum <= eps*errden) Then
Exit restart
End If
If (itrace_ > 0) then
if (((itx==0).or.(mod(itx,itrace_)==0)).and.(me == 0)) &
& write(*,'(a,i4,3(2x,es10.4))') 'bicgstab: ',&
& itx,errnum,eps*errden
Endif
if (errnum <= eps*errden) exit restart
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,1,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(q,r,desc_a,info)
If (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' RHO:',rho
If (rho==dzero) Then
If (debug_level >= psb_debug_ext_) &
if (rho==dzero) then
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Iteration breakdown R',rho
Exit iteration
Endif
exit iteration
endif
If (it==1) Then
Call psb_geaxpby(done,r,dzero,p,desc_a,info)
Else
if (it==1) then
call psb_geaxpby(done,r,dzero,p,desc_a,info)
else
beta = (rho/rho_old)*(alpha/omega)
Call psb_geaxpby(-omega,v,done,p,desc_a,info)
Call psb_geaxpby(done,r,beta,p,desc_a,info)
call psb_geaxpby(-omega,v,done,p,desc_a,info)
call psb_geaxpby(done,r,beta,p,desc_a,info)
End If
#ifdef MPE_KRYLOV
imerr = MPE_Log_event( ifctb, 0, "st PREC" )
#endif
Call psb_precaply(prec,p,f,desc_a,info,work=aux)
call psb_precaply(prec,p,f,desc_a,info,work=aux)
#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,f,dzero,v,desc_a,info,&
call psb_spmm(done,a,f,dzero,v,desc_a,info,&
& work=aux)
#ifdef MPE_KRYLOV
imerr = MPE_Log_event( imme, 0, "ed SPMM" )
#endif
sigma = psb_gedot(q,v,desc_a,info)
If (sigma==dzero) Then
If (debug_level >= psb_debug_ext_) &
if (sigma==dzero) 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
exit iteration
endif
alpha = rho/sigma
Call psb_geaxpby(done,r,dzero,s,desc_a,info)
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)
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
@@ -388,7 +376,7 @@ Subroutine psb_dcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
#ifdef MPE_KRYLOV
imerr = MPE_Log_event( ifctb, 0, "st PREC" )
#endif
Call psb_precaply(prec,s,z,desc_a,info,work=aux)
call psb_precaply(prec,s,z,desc_a,info,work=aux)
#ifdef MPE_KRYLOV
imerr = MPE_Log_event( ifcte, 0, "ed PREC" )
#endif
@@ -412,91 +400,79 @@ Subroutine psb_dcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
end if
sigma = psb_gedot(t,t,desc_a,info)
If (sigma==dzero) Then
If (debug_level >= psb_debug_ext_) &
if (sigma==dzero) then
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Iteration breakdown S2', sigma
Exit iteration
Endif
exit iteration
endif
tau = psb_gedot(t,s,desc_a,info)
omega = tau/sigma
If (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Omega:',omega
If (omega==dzero) Then
If (debug_level >= psb_debug_ext_) &
if (omega==dzero) then
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Iteration breakdown O',omega
Exit iteration
Endif
exit iteration
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)
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 (istop_ == 1) Then
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
else if (istop_ == 2) then
rni = psb_genrm2(r,desc_a,info)
errnum = rni
errden = bn2
Endif
If (errnum <= eps*errden) Then
Exit restart
End If
If (itx.Ge.litmax) Exit restart
endif
If (itrace_ > 0) then
if ((mod(itx,itrace_)==0).and.(me == 0)) &
& write(*,'(a,i4,3(2x,es10.4))') &
& 'bicgstab: ',itx,errnum,eps*errden
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)
End Do iteration
End Do restart
end do iteration
end do restart
If (itrace_ > 0) then
if (me == 0) write(*,'(a,i4,3(2x,es10.4))') 'bicgstab: ',&
& itx,errnum,eps*errden
Endif
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,1,errnum,errden,eps)
If (Present(err)) then
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).and.(me==0)) Then
write(debug_unit,*) 'BI-CGSTAB failed to converge to ',eps*errden,&
& ' in ',ITX,' iterations. '
End If
end if
if (present(iter)) iter = itx
if (errnum > eps*errden) &
& call end_log(methdname,me,itx,errnum,errden,eps)
Deallocate(aux)
Call psb_gefree(wwrk,desc_a,info)
! restore external global coherence behaviour
call psb_restore_coher(ictxt,isvch)
#ifdef MPE_KRYLOV
imerr = MPE_Log_event( istpe, 0, "ed CGSTAB" )
#endif
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
end if
#ifdef MPE_KRYLOV
imerr = MPE_Log_event( istpe, 0, "ed CGSTAB" )
#endif
! restore external global coherence behaviour
call psb_restore_coher(ictxt,isvch)
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
if (err_act == psb_act_abort_) then
call psb_error(ictxt)
return
end if
krylov/psb_dcgstabl.f90
+106 -127
View File
@@ -105,6 +105,7 @@
Subroutine psb_dcgstabl(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,istop)
use psb_base_mod
use psb_prec_mod
use psb_krylov_mod, psb_protect_name => psb_dcgstabl
implicit none
!!$ parameters
@@ -133,6 +134,7 @@ Subroutine psb_dcgstabl(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,is
& omega
real(kind(1.d0)) :: errnum, errden
character(len=20) :: name
character(len=*), parameter :: methdname='BiCGStab(L)'
info = 0
name = 'psb_dcgstabl'
@@ -168,29 +170,29 @@ Subroutine psb_dcgstabl(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,is
goto 9999
endif
If (Present(itmax)) Then
if (present(itmax)) then
litmax = itmax
Else
else
litmax = 1000
Endif
endif
If (Present(itrace)) Then
if (present(itrace)) then
itrace_ = itrace
Else
else
itrace_ = 0
End If
end if
If (Present(irst)) Then
if (present(irst)) then
nl = irst
If (debug_level >= psb_debug_ext_) &
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'present: irst: ',irst,nl
Else
else
nl = 1
If (debug_level >= psb_debug_ext_) &
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' not present: irst: ',irst,nl
Endif
endif
if (nl <=0 ) then
info=5001
int_err(1)=nl
@@ -213,25 +215,25 @@ Subroutine psb_dcgstabl(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,is
end if
naux=4*n_col
Allocate(aux(naux),gamma(0:nl),gamma1(nl),&
allocate(aux(naux),gamma(0:nl),gamma1(nl),&
&gamma2(nl),taum(nl,nl),sigma(nl), stat=info)
If (info.Ne.0) Then
if (info /= 0) then
info=4000
call psb_errpush(info,name)
goto 9999
End If
end if
if (info == 0) Call psb_geall(wwrk,desc_a,info,n=10)
if (info == 0) Call psb_geall(uh,desc_a,info,n=nl+1)
if (info == 0) Call psb_geall(rh,desc_a,info,n=nl+1)
if (info == 0) Call psb_geasb(wwrk,desc_a,info)
if (info == 0) Call psb_geasb(uh,desc_a,info)
if (info == 0) Call psb_geasb(rh,desc_a,info)
if (info.ne.0) Then
if (info /= 0) then
info=4011
call psb_errpush(info,name)
goto 9999
End If
end if
q => wwrk(:,1)
r => wwrk(:,2)
@@ -255,40 +257,40 @@ Subroutine psb_dcgstabl(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,is
endif
errnum = dzero
errden = done
if (info.ne.0) Then
if (info /= 0) then
info=4011
call psb_errpush(info,name)
goto 9999
End If
end if
itx = 0
restart: Do
restart: do
!!$
!!$ r0 = b-ax0
!!$
If (debug_level >= psb_debug_ext_) &
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),' restart: ',itx,it
If (itx.Ge.litmax) Exit restart
if (itx >= litmax) 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,b,dzero,r,desc_a,info)
call psb_spmm(-done,a,x,done,r,desc_a,info,work=aux)
call psb_precaply(prec,r,desc_a,info)
Call psb_geaxpby(done,r,dzero,rt0,desc_a,info)
Call psb_geaxpby(done,r,dzero,rh(:,0),desc_a,info)
Call psb_geaxpby(dzero,r,dzero,uh(:,0),desc_a,info)
if (info.ne.0) Then
call psb_geaxpby(done,r,dzero,rt0,desc_a,info)
call psb_geaxpby(done,r,dzero,rh(:,0),desc_a,info)
call psb_geaxpby(dzero,r,dzero,uh(:,0),desc_a,info)
if (info /= 0) then
info=4011
call psb_errpush(info,name)
goto 9999
End If
end if
rho = done
alpha = dzero
omega = done
If (debug_level >= psb_debug_ext_) &
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' on entry to amax: b: ',Size(b)
@@ -302,122 +304,107 @@ Subroutine psb_dcgstabl(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,is
errnum = rni
errden = bn2
endif
if (info.ne.0) Then
if (info /= 0) Then
info=4011
call psb_errpush(info,name)
goto 9999
End If
end if
If (errnum <= eps*errden) Then
Exit restart
End If
If (itrace_ > 0) then
if ((mod(itx,itrace_)==0).and.(me == 0))&
& write(*,'(a,i4,3(2x,es10.4))') 'bicgstab(l): ',itx,errnum,eps*errden
end If
if (errnum <= eps*errden) exit restart
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,eps)
iteration: Do
iteration: do
it = it + nl
itx = itx + nl
rho = -omega*rho
If (debug_level >= psb_debug_ext_) &
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' iteration: ',itx, rho,rh(1,0)
Do j = 0, nl -1
do j = 0, nl -1
If (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),'bicg part: ',j, nl
rho_old = rho
rho = psb_gedot(rh(:,j),rt0,desc_a,info)
If (rho==dzero) Then
If (debug_level >= psb_debug_ext_) &
if (rho==dzero) then
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' bi-cgstab iteration breakdown r',rho
Exit iteration
Endif
exit iteration
endif
beta = alpha*rho/rho_old
If (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' bicg part: ',alpha,beta,rho,rho_old
rho_old = rho
Call psb_geaxpby(done,rh(:,0:j),-beta,uh(:,0:j),desc_a,info)
If (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' bicg part: ',rh(1,0),beta
Call psb_spmm(done,a,uh(:,j),dzero,uh(:,j+1),desc_a,info,work=aux)
call psb_geaxpby(done,rh(:,0:j),-beta,uh(:,0:j),desc_a,info)
call psb_spmm(done,a,uh(:,j),dzero,uh(:,j+1),desc_a,info,work=aux)
call psb_precaply(prec,uh(:,j+1),desc_a,info)
gamma(j) = psb_gedot(uh(:,j+1),rt0,desc_a,info)
If (gamma(j)==dzero) Then
If (debug_level >= psb_debug_ext_) &
if (gamma(j)==dzero) then
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' bi-cgstab iteration breakdown s2',gamma(j)
Exit iteration
Endif
exit iteration
endif
alpha = rho/gamma(j)
If (debug_level >= psb_debug_ext_) &
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' bicg part: alpha=r/g ',alpha,rho,gamma(j)
Call psb_geaxpby(-alpha,uh(:,1:j+1),done,rh(:,0:j),desc_a,info)
Call psb_geaxpby(alpha,uh(:,0),done,x,desc_a,info)
Call psb_spmm(done,a,rh(:,j),dzero,rh(:,j+1),desc_a,info,work=aux)
call psb_geaxpby(-alpha,uh(:,1:j+1),done,rh(:,0:j),desc_a,info)
call psb_geaxpby(alpha,uh(:,0),done,x,desc_a,info)
call psb_spmm(done,a,rh(:,j),dzero,rh(:,j+1),desc_a,info,work=aux)
call psb_precaply(prec,rh(:,j+1),desc_a,info)
Enddo
enddo
Do j=1, nl
If (debug_level >= psb_debug_ext_) &
do j=1, nl
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' mod g-s part: ',j, nl,rh(1,0)
Do i=1, j-1
do i=1, j-1
taum(i,j) = psb_gedot(rh(:,i),rh(:,j),desc_a,info)
taum(i,j) = taum(i,j)/sigma(i)
Call psb_geaxpby(-taum(i,j),rh(:,i),done,rh(:,j),desc_a,info)
Enddo
call psb_geaxpby(-taum(i,j),rh(:,i),done,rh(:,j),desc_a,info)
enddo
sigma(j) = psb_gedot(rh(:,j),rh(:,j),desc_a,info)
gamma1(j) = psb_gedot(rh(:,0),rh(:,j),desc_a,info)
If (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' mod g-s part: gamma1 ',gamma1(j), sigma(j)
gamma1(j) = gamma1(j)/sigma(j)
Enddo
enddo
gamma(nl) = gamma1(nl)
omega = gamma(nl)
Do j=nl-1,1,-1
do j=nl-1,1,-1
gamma(j) = gamma1(j)
Do i=j+1,nl
do i=j+1,nl
gamma(j) = gamma(j) - taum(j,i) * gamma(i)
Enddo
Enddo
If (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' first solve: ', gamma(:)
Do j=1,nl-1
gamma2(j) = gamma(j+1)
Do i=j+1,nl-1
gamma2(j) = gamma2(j) + taum(j,i) * gamma(i+1)
Enddo
Enddo
If (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' second solve: ', gamma(:)
Call psb_geaxpby(gamma(1),rh(:,0),done,x,desc_a,info)
Call psb_geaxpby(-gamma1(nl),rh(:,nl),done,rh(:,0),desc_a,info)
Call psb_geaxpby(-gamma(nl),uh(:,nl),done,uh(:,0),desc_a,info)
enddo
enddo
Do j=1, nl-1
Call psb_geaxpby(-gamma(j),uh(:,j),done,uh(:,0),desc_a,info)
Call psb_geaxpby(gamma2(j),rh(:,j),done,x,desc_a,info)
Call psb_geaxpby(-gamma1(j),rh(:,j),done,rh(:,0),desc_a,info)
Enddo
do j=1,nl-1
gamma2(j) = gamma(j+1)
do i=j+1,nl-1
gamma2(j) = gamma2(j) + taum(j,i) * gamma(i+1)
enddo
enddo
call psb_geaxpby(gamma(1),rh(:,0),done,x,desc_a,info)
call psb_geaxpby(-gamma1(nl),rh(:,nl),done,rh(:,0),desc_a,info)
call psb_geaxpby(-gamma(nl),uh(:,nl),done,uh(:,0),desc_a,info)
do j=1, nl-1
call psb_geaxpby(-gamma(j),uh(:,j),done,uh(:,0),desc_a,info)
call psb_geaxpby(gamma2(j),rh(:,j),done,x,desc_a,info)
call psb_geaxpby(-gamma1(j),rh(:,j),done,rh(:,0),desc_a,info)
enddo
if (istop_ == 1) then
rni = psb_geamax(rh(:,0),desc_a,info)
@@ -430,54 +417,46 @@ Subroutine psb_dcgstabl(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,is
errden = bn2
endif
If (errnum <= eps*errden) Then
Exit restart
End If
If (itx.Ge.litmax) Exit restart
If (itrace_ > 0) then
if ((mod(itx,itrace_)==0).and.(me == 0))&
& write(*,'(a,i4,3(2x,es10.4))') 'bicgstab(l): ',itx,errnum,eps*errden
end If
if (errnum <= eps*errden) exit restart
if (itx >= litmax) exit restart
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,eps)
End Do iteration
End Do restart
end do iteration
end do restart
If (itrace_ > 0) then
if (me == 0) write(*,'(a,i4,3(2x,es10.4))') 'bicgstab(l): ',itx,errnum,eps*errden
end If
If (Present(err)) then
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).and.(me==0)) Then
write(debug_unit,*) 'bi-cgstabl failed to converge to ',eps*errden,&
& ' in ',itx,' iterations '
End If
end if
Deallocate(aux)
Call psb_gefree(wwrk,desc_a,info)
Call psb_gefree(uh,desc_a,info)
Call psb_gefree(rh,desc_a,info)
! restore external global coherence behaviour
call psb_restore_coher(ictxt,isvch)
if (present(iter)) iter = itx
if (errnum > eps*errden) &
& call end_log(methdname,me,itx,errnum,errden,eps)
deallocate(aux,stat=info)
if (info == 0) call psb_gefree(wwrk,desc_a,info)
if (info == 0) call psb_gefree(uh,desc_a,info)
if (info == 0) call psb_gefree(rh,desc_a,info)
if(info/=0) then
call psb_errpush(info,name)
goto 9999
end if
! restore external global coherence behaviour
call psb_restore_coher(ictxt,isvch)
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
krylov/psb_drgmres.f90
+68 -79
View File
@@ -105,9 +105,10 @@
! irst - integer(optional) Input: restart parameter
!
!
Subroutine psb_drgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,istop)
subroutine psb_drgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,istop)
use psb_base_mod
use psb_prec_mod
use psb_krylov_mod, psb_protect_name => psb_drgmres
implicit none
!!$ Parameters
@@ -134,7 +135,8 @@ Subroutine psb_drgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist
Real(Kind(1.d0)) :: rni, xni, bni, ani,bn2, dt
real(kind(1.d0)), external :: dnrm2
real(kind(1.d0)) :: errnum, errden
character(len=20) :: name
character(len=20) :: name
character(len=*), parameter :: methdname='RGMRES'
info = 0
name = 'psb_dgmres'
@@ -169,29 +171,29 @@ Subroutine psb_drgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist
goto 9999
endif
If (Present(itmax)) Then
if (present(itmax)) then
litmax = itmax
Else
else
litmax = 1000
Endif
endif
If (Present(itrace)) Then
if (present(itrace)) then
itrace_ = itrace
Else
else
itrace_ = 0
End If
end if
If (Present(irst)) Then
if (present(irst)) then
nl = irst
If (debug_level >= psb_debug_ext_) &
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' present: irst: ',irst,nl
Else
else
nl = 10
If (debug_level >= psb_debug_ext_) &
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' not present: irst: ',irst,nl
Endif
endif
if (nl <=0 ) then
info=5001
int_err(1)=nl
@@ -215,22 +217,22 @@ Subroutine psb_drgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist
naux=4*n_col
Allocate(aux(naux),h(nl+1,nl+1),&
allocate(aux(naux),h(nl+1,nl+1),&
&c(nl+1),s(nl+1),rs(nl+1), rst(nl+1),stat=info)
if (info == 0) Call psb_geall(v,desc_a,info,n=nl+1)
if (info == 0) Call psb_geall(w,desc_a,info)
if (info == 0) Call psb_geall(w1,desc_a,info)
if (info == 0) Call psb_geall(xt,desc_a,info)
if (info == 0) Call psb_geasb(v,desc_a,info)
if (info == 0) Call psb_geasb(w,desc_a,info)
if (info == 0) Call psb_geasb(w1,desc_a,info)
if (info == 0) Call psb_geasb(xt,desc_a,info)
if (info.ne.0) Then
if (info == 0) call psb_geall(v,desc_a,info,n=nl+1)
if (info == 0) call psb_geall(w,desc_a,info)
if (info == 0) call psb_geall(w1,desc_a,info)
if (info == 0) call psb_geall(xt,desc_a,info)
if (info == 0) call psb_geasb(v,desc_a,info)
if (info == 0) call psb_geasb(w,desc_a,info)
if (info == 0) call psb_geasb(w1,desc_a,info)
if (info == 0) call psb_geasb(xt,desc_a,info)
if (info /= 0) then
info=4011
call psb_errpush(info,name)
goto 9999
End If
end if
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Size of V,W,W1 ',size(v),size(v,1),&
@@ -247,47 +249,47 @@ Subroutine psb_drgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist
endif
errnum = dzero
errden = done
if (info.ne.0) Then
if (info /= 0) then
info=4011
call psb_errpush(info,name)
goto 9999
End If
end if
itx = 0
restart: Do
restart: do
! compute r0 = b-ax0
! check convergence
! compute v1 = r0/||r0||_2
If (debug_level >= psb_debug_ext_) &
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' restart: ',itx,it
it = 0
Call psb_geaxpby(done,b,dzero,v(:,1),desc_a,info)
if (info.ne.0) Then
call psb_geaxpby(done,b,dzero,v(:,1),desc_a,info)
if (info /= 0) then
info=4011
call psb_errpush(info,name)
goto 9999
End If
end if
Call psb_spmm(-done,a,x,done,v(:,1),desc_a,info,work=aux)
if (info.ne.0) Then
call psb_spmm(-done,a,x,done,v(:,1),desc_a,info,work=aux)
if (info /= 0) then
info=4011
call psb_errpush(info,name)
goto 9999
End If
end if
rs(1) = psb_genrm2(v(:,1),desc_a,info)
rs(2:) = dzero
if (info.ne.0) Then
if (info /= 0) then
info=4011
call psb_errpush(info,name)
goto 9999
End If
end if
scal=done/rs(1) ! rs(1) MIGHT BE VERY SMALL - USE DSCAL TO DEAL WITH IT?
If (debug_level >= psb_debug_ext_) &
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' on entry to amax: b: ',Size(b),rs(1),scal
@@ -304,24 +306,20 @@ Subroutine psb_drgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist
errnum = rni
errden = bn2
endif
if (info.ne.0) Then
if (info /= 0) then
info=4011
call psb_errpush(info,name)
goto 9999
End If
end if
If (errnum <= eps*errden) Then
Exit restart
End If
if (errnum <= eps*errden) exit restart
If (itrace_ > 0) then
if ((mod(itx,itrace_)==0).and.(me == 0))&
& write(*,'(a,i4,3(2x,es10.4))') 'gmres(l): ',itx,errnum,eps*errden
end If
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,eps)
v(:,1) = v(:,1) * scal
If (itx.Ge.litmax) Exit restart
if (itx >= litmax) exit restart
!
! inner iterations
@@ -331,7 +329,7 @@ Subroutine psb_drgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist
itx = itx + 1
call psb_precaply(prec,v(:,i),w1,desc_a,info)
Call psb_spmm(done,a,w1,dzero,w,desc_a,info,work=aux)
call psb_spmm(done,a,w1,dzero,w,desc_a,info,work=aux)
!
do k = 1, i
@@ -403,7 +401,7 @@ Subroutine psb_drgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist
errden = bn2
endif
If (errnum <= eps*errden) Then
if (errnum <= eps*errden) then
if (istop_ == 1) then
x = xt
@@ -427,12 +425,10 @@ Subroutine psb_drgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist
end if
If (itrace_ > 0) then
if ((mod(itx,itrace_)==0).and.(me == 0))&
& write(*,'(a,i4,3(2x,es10.4))') 'gmres(l): ',itx,errnum,eps*errden
end If
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,eps)
end Do inner
end do inner
if (istop_ == 1) then
x = xt
@@ -452,53 +448,47 @@ Subroutine psb_drgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist
call psb_geaxpby(done,w,done,x,desc_a,info)
end if
End Do restart
If (itrace_ > 0) then
if (me == 0) write(*,'(a,i4,3(2x,es10.4))') 'gmres(l): ',itx,errnum,eps*errden
end If
end do restart
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,1,errnum,errden,eps)
If (Present(err)) then
if (present(err)) then
if (errden /= dzero) then
err = errnum/errden
else
err = errnum
end if
end If
end if
If (Present(iter)) iter = itx
If ((errnum > eps*errden).and.(me==0)) Then
write(debug_unit,*) 'gmresr(l) failed to converge to ',eps*errden,&
& ' in ',itx,' iterations '
End If
Deallocate(aux,h,c,s,rs,rst, stat=info)
Call psb_gefree(v,desc_a,info)
Call psb_gefree(w,desc_a,info)
Call psb_gefree(w1,desc_a,info)
Call psb_gefree(xt,desc_a,info)
! restore external global coherence behaviour
call psb_restore_coher(ictxt,isvch)
if (present(iter)) iter = itx
if (errnum > eps*errden) &
& call end_log(methdname,me,itx,errnum,errden,eps)
deallocate(aux,h,c,s,rs,rst, stat=info)
if (info == 0) call psb_gefree(v,desc_a,info)
if (info == 0) call psb_gefree(w,desc_a,info)
if (info == 0) call psb_gefree(w1,desc_a,info)
if (info == 0) call psb_gefree(xt,desc_a,info)
if (info /= 0) then
info=4011
call psb_errpush(info,name)
goto 9999
end if
! restore external global coherence behaviour
call psb_restore_coher(ictxt,isvch)
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
contains
function safe_dn2(a,b)
real(kind(1.d0)), intent(in) :: a, b
@@ -513,8 +503,7 @@ contains
endif
return
end function safe_dn2
End Subroutine psb_drgmres
end subroutine psb_drgmres
krylov/psb_krylov_mod.f90
+26
View File
@@ -423,6 +423,32 @@ contains
end subroutine psb_zkrylov
subroutine log_conv(methdname,me,itx,itrace,errnum,errden,eps)
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
Endif
end subroutine log_conv
subroutine end_log(methdname,me,itx,errnum,errden,eps)
character(len=*), intent(in) :: methdname
integer, intent(in) :: me, itx
real(kind(1.d0)), intent(in) :: 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))'
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
end if
end subroutine end_log
end module psb_krylov_mod
krylov/psb_zcgs.f90
+58 -73
View File
@@ -94,6 +94,7 @@
Subroutine psb_zcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
use psb_base_mod
use psb_prec_mod
use psb_krylov_mod, psb_protect_name => psb_zcgs
implicit none
!!$ parameters
@@ -121,6 +122,7 @@ Subroutine psb_zcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
complex(Kind(1.d0)) :: alpha, beta, rho, rho_old, sigma
real(kind(1.d0)) :: errnum, errden
character(len=20) :: name
character(len=*), parameter :: methdname='CGS'
info = 0
name = 'psb_zcgs'
@@ -197,11 +199,11 @@ Subroutine psb_zcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
litmax = 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)
@@ -226,7 +228,7 @@ Subroutine psb_zcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
!!$
!!$ r0 = b-ax0
!!$
If (itx.Ge.litmax) Exit restart
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)
@@ -258,68 +260,63 @@ Subroutine psb_zcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
goto 9999
end if
If (errnum <= eps*errden) Then
Exit restart
End If
If (itrace_ > 0) then
if ((mod(itx,itrace_)==0).and.(me == 0))&
& write(*,'(a,i4,3(2x,es10.4))') 'cgs: ',itx,errnum,eps*errden
end If
if (errnum <= eps*errden) exit restart
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==zzero) Then
If (debug_level >= psb_debug_ext_) &
if (rho==zzero) then
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' iteration breakdown r',rho
Exit iteration
Endif
exit iteration
endif
If (it==1) Then
Call psb_geaxpby(zone,r,zzero,uv,desc_a,info)
Call psb_geaxpby(zone,r,zzero,p,desc_a,info)
Else
if (it==1) then
call psb_geaxpby(zone,r,zzero,uv,desc_a,info)
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)
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)
end if
End If
call psb_precaply(prec,p,f,desc_a,info,work=aux)
Call psb_precaply(prec,p,f,desc_a,info,work=aux)
Call psb_spmm(zone,a,f,zzero,v,desc_a,info,&
call psb_spmm(zone,a,f,zzero,v,desc_a,info,&
& work=aux)
sigma = psb_gedot(rt,v,desc_a,info)
If (sigma==zzero) Then
If (debug_level >= psb_debug_ext_) &
if (sigma==zzero) then
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' iteration breakdown s1', sigma
Exit iteration
Endif
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_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)
call psb_precaply(prec,s,z,desc_a,info,work=aux)
Call psb_geaxpby(alpha,z,zone,x,desc_a,info)
call psb_geaxpby(alpha,z,zone,x,desc_a,info)
Call psb_spmm(zone,a,z,zzero,qt,desc_a,info,&
call psb_spmm(zone,a,z,zzero,qt,desc_a,info,&
& work=aux)
Call psb_geaxpby(-alpha,qt,zone,r,desc_a,info)
call psb_geaxpby(-alpha,qt,zone,r,desc_a,info)
if (istop_ == 1) then
@@ -333,59 +330,47 @@ Subroutine psb_zcgs(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
errden = bn2
endif
If (errnum <= eps*errden) Then
Exit restart
End If
If (itx.Ge.litmax) Exit restart
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 (itrace_ > 0) then
if ((mod(itx,itrace_)==0).and.(me == 0))&
& write(*,'(a,i4,3(2x,es10.4))') 'cgs: ',itx,errnum,eps*errden
end If
end do iteration
end do restart
End Do iteration
End Do restart
If (itrace_ > 0) then
if ((mod(itx,itrace_)==0).and.(me == 0))&
& write(*,'(a,i4,3(2x,es10.4))') 'cgs: ',itx,errnum,eps*errden
end If
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,1,errnum,errden,eps)
If (Present(err)) then
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).and.(me==0)) Then
write(debug_unit,*) 'cgs failed to converge to ',eps*errden,&
& ' in ',itx,' iterations '
End If
end if
if (present(iter)) iter = itx
if (errnum > eps*errden) &
& call end_log(methdname,me,itx,errnum,errden,eps)
Deallocate(aux)
Call psb_gefree(wwrk,desc_a,info)
! restore external global coherence behaviour
call psb_restore_coher(ictxt,isvch)
if(info/=0) then
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
end if
! restore external global coherence behaviour
call psb_restore_coher(ictxt,isvch)
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
End Subroutine psb_zcgs
end subroutine psb_zcgs
krylov/psb_zcgstab.f90
+92 -115
View File
@@ -92,9 +92,10 @@
! where r is the (preconditioned, recursive
! estimate of) residual.
!
Subroutine psb_zcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
subroutine psb_zcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
use psb_base_mod
use psb_prec_mod
use psb_krylov_mod, psb_protect_name => psb_zcgstab
Implicit None
!!$ parameters
Type(psb_zspmat_type), Intent(in) :: a
@@ -123,7 +124,8 @@ Subroutine psb_zcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
real(kind(1.d0)) :: errnum, errden
!!$ 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
character(len=20) :: name
character(len=*), parameter :: methdname='BiCGStab'
info = 0
name = 'psb_zcgstab'
@@ -189,17 +191,17 @@ Subroutine psb_zcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
T => WWRK(:,7)
Z => WWRK(:,8)
If (Present(itmax)) Then
if (present(itmax)) then
litmax = itmax
Else
else
litmax = 1000
Endif
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)
@@ -222,19 +224,19 @@ Subroutine psb_zcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
!!$
!!$ r0 = b-Ax0
!!$
If (itx >= litmax) Exit restart
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) Then
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) then
info=4011
call psb_errpush(info,name)
goto 9999
End If
end if
rho = zzero
If (debug_level >= psb_debug_ext_) &
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' On entry to AMAX: B: ',Size(b)
@@ -242,52 +244,42 @@ Subroutine psb_zcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
! Must always provide norm of R into RNI below for first check on
! residual
!
If (istop_ == 1) Then
if (istop_ == 1) then
rni = psb_geamax(r,desc_a,info)
xni = psb_geamax(x,desc_a,info)
Else If (istop_ == 2) Then
else if (istop_ == 2) then
rni = psb_genrm2(r,desc_a,info)
Endif
endif
errnum = dzero
errden = done
if (info /= 0) Then
if (info /= 0) then
info=4011
call psb_errpush(info,name)
goto 9999
End If
end if
If (itx == 0) Then
if (itx == 0) then
rn0 = rni
End If
If (rn0 == 0.d0 ) Then
If (itrace_ > 0 ) Then
If (me == 0) Write(*,*) 'BiCGSTAB: ',itx,rn0
Endif
Exit restart
End If
end if
if (rn0 == 0.d0 ) exit restart
If (istop_ == 1) Then
if (istop_ == 1) then
xni = psb_geamax(x,desc_a,info)
errnum = rni
errden = (ani*xni+bni)
Else If (istop_ == 2) Then
else if (istop_ == 2) then
errnum = rni
errden = bn2
Endif
if (info /= 0) Then
endif
if (info /= 0) then
info=4011
call psb_errpush(info,name)
goto 9999
End If
end if
If (errnum <= eps*errden) Then
Exit restart
End If
If (itrace_ > 0) then
if (((itx==0).or.(mod(itx,itrace_)==0)).and.(me == 0)) &
& write(*,'(a,i4,3(2x,es10.4))') 'bicgstab: ',itx,errnum,eps*errden
end If
if (errnum <= eps*errden) exit restart
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,eps)
iteration: Do
it = it + 1
@@ -298,148 +290,133 @@ Subroutine psb_zcgstab(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,istop)
rho_old = rho
rho = psb_gedot(q,r,desc_a,info)
If (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' RHO:',rho
If (rho==zzero) Then
If (debug_level >= psb_debug_ext_) &
if (rho==zzero) then
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Iteration breakdown R',rho
Exit iteration
Endif
exit iteration
endif
If (it==1) Then
Call psb_geaxpby(zone,r,zzero,p,desc_a,info)
Else
if (it==1) then
call psb_geaxpby(zone,r,zzero,p,desc_a,info)
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)
End If
call psb_geaxpby(-omega,v,zone,p,desc_a,info)
call psb_geaxpby(zone,r,beta,p,desc_a,info)
end if
Call psb_precaply(prec,p,f,desc_a,info,work=aux)
call psb_precaply(prec,p,f,desc_a,info,work=aux)
Call psb_spmm(zone,a,f,zzero,v,desc_a,info,&
call psb_spmm(zone,a,f,zzero,v,desc_a,info,&
& work=aux)
sigma = psb_gedot(q,v,desc_a,info)
If (sigma==zzero) Then
If (debug_level >= psb_debug_ext_) &
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_) &
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_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_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_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_spmm(zone,a,z,zzero,t,desc_a,info,&
& work=aux)
if(info /= 0) then
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_spmm')
goto 9999
end if
sigma = psb_gedot(t,t,desc_a,info)
If (sigma==zzero) Then
If (debug_level >= psb_debug_ext_) &
if (sigma==zzero) then
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Iteration breakdown S2', sigma
Exit iteration
Endif
exit iteration
endif
tau = psb_gedot(t,s,desc_a,info)
omega = tau/sigma
If (omega==zzero) Then
If (debug_level >= psb_debug_ext_) &
if (omega==zzero) then
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Iteration breakdown O',omega
Exit iteration
Endif
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)
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
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
else if (istop_ == 2) then
rni = psb_genrm2(r,desc_a,info)
errnum = rni
errden = bn2
Endif
endif
If (errnum <= eps*errden) Then
Exit restart
End If
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 (itx.Ge.litmax) Exit restart
end do iteration
end do restart
If (itrace_ > 0) then
if ((mod(itx,itrace_)==0).and.(me == 0)) &
& write(*,'(a,i4,3(2x,es10.4))') &
& 'bicgstab: ',itx,errnum,eps*errden
Endif
End Do iteration
End Do restart
If (itrace_ > 0) then
if (me == 0) write(*,'(a,i4,3(2x,es10.4))') 'bicgstab: ',itx,errnum,eps*errden
Endif
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,1,errnum,errden,eps)
If (Present(err)) then
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).and.(me==0)) Then
write(debug_unit,*) 'BI-CGSTAB failed to converge to ',eps*errden,&
& ' in ',ITX,' iterations. '
End If
end if
if (present(iter)) iter = itx
if (errnum > eps*errden) &
& call end_log(methdname,me,itx,errnum,errden,eps)
Deallocate(aux)
Call psb_gefree(wwrk,desc_a,info)
! restore external global coherence behaviour
call psb_restore_coher(ictxt,isvch)
!!$ imerr = MPE_Log_event( istpe, 0, "ed CGSTAB" )
if(info/=0) then
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
end if
! restore external global coherence behaviour
call psb_restore_coher(ictxt,isvch)
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
if (err_act == psb_act_abort_) then
call psb_error(ictxt)
return
end if
krylov/psb_zrgmres.f90
+55 -64
View File
@@ -107,6 +107,7 @@
Subroutine psb_zrgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,istop)
use psb_base_mod
use psb_prec_mod
use psb_krylov_mod, psb_protect_name => psb_zrgmres
implicit none
!!$ Parameters
@@ -134,7 +135,8 @@ Subroutine psb_zrgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist
Real(Kind(1.d0)) :: rni, xni, bni, ani,bn2
real(kind(1.d0)), external :: dznrm2
real(kind(1.d0)) :: errnum, errden
character(len=20) :: name
character(len=20) :: name
character(len=*), parameter :: methdname='RGMRES'
info = 0
name = 'psb_zgmres'
@@ -169,29 +171,29 @@ Subroutine psb_zrgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist
goto 9999
endif
If (Present(itmax)) Then
if (present(itmax)) then
litmax = itmax
Else
else
litmax = 1000
Endif
endif
If (Present(itrace)) Then
if (present(itrace)) then
itrace_ = itrace
Else
else
itrace_ = 0
End If
end if
If (Present(irst)) Then
if (present(irst)) then
nl = irst
If (debug_level >= psb_debug_ext_) &
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' present: irst: ',irst,nl
Else
else
nl = 10
If (debug_level >= psb_debug_ext_) &
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' not present: irst: ',irst,nl
Endif
endif
if (nl <=0 ) then
info=5001
int_err(1)=nl
@@ -215,7 +217,7 @@ Subroutine psb_zrgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist
naux=4*n_col
Allocate(aux(naux),h(nl+1,nl+1),&
allocate(aux(naux),h(nl+1,nl+1),&
&c(nl+1),s(nl+1),rs(nl+1), rst(nl+1),stat=info)
if (info == 0) Call psb_geall(v,desc_a,info,n=nl+1)
@@ -226,11 +228,11 @@ Subroutine psb_zrgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist
if (info == 0) Call psb_geasb(w,desc_a,info)
if (info == 0) Call psb_geasb(w1,desc_a,info)
if (info == 0) Call psb_geasb(xt,desc_a,info)
if (info.ne.0) Then
if (info /= 0) then
info=4011
call psb_errpush(info,name)
goto 9999
End If
end if
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Size of V,W,W1 ',size(v),size(v,1),&
@@ -247,11 +249,11 @@ Subroutine psb_zrgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist
endif
errnum = dzero
errden = done
if (info.ne.0) Then
if (info /= 0) then
info=4011
call psb_errpush(info,name)
goto 9999
End If
end if
itx = 0
restart: Do
@@ -260,19 +262,19 @@ Subroutine psb_zrgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist
! check convergence
! compute v1 = r0/||r0||_2
If (debug_level >= psb_debug_ext_) &
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' restart: ',itx,it
it = 0
Call psb_geaxpby(zone,b,zzero,v(:,1),desc_a,info)
if (info.ne.0) Then
call psb_geaxpby(zone,b,zzero,v(:,1),desc_a,info)
if (info /= 0) then
info=4011
call psb_errpush(info,name)
goto 9999
End If
end if
Call psb_spmm(-zone,a,x,zone,v(:,1),desc_a,info,work=aux)
if (info.ne.0) Then
call psb_spmm(-zone,a,x,zone,v(:,1),desc_a,info,work=aux)
if (info /= 0) Then
info=4011
call psb_errpush(info,name)
goto 9999
@@ -280,14 +282,14 @@ Subroutine psb_zrgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist
rs(1) = psb_genrm2(v(:,1),desc_a,info)
rs(2:) = zzero
if (info.ne.0) Then
if (info /= 0) Then
info=4011
call psb_errpush(info,name)
goto 9999
End If
end if
scal=done/rs(1) ! rs(1) MIGHT BE VERY SMALL - USE DSCAL TO DEAL WITH IT?
If (debug_level >= psb_debug_ext_) &
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' on entry to amax: b: ',Size(b),rs(1),scal
@@ -304,24 +306,20 @@ Subroutine psb_zrgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist
errnum = rni
errden = bn2
endif
if (info.ne.0) Then
if (info /= 0) then
info=4011
call psb_errpush(info,name)
goto 9999
End If
end if
If (errnum <= eps*errden) Then
Exit restart
End If
if (errnum <= eps*errden) exit restart
If (itrace_ > 0) then
if ((mod(itx,itrace_)==0).and.(me == 0))&
& write(*,'(a,i4,3(2x,es10.4))') 'gmres(l): ',itx,errnum,eps*errden
end If
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,eps)
v(:,1) = v(:,1) * scal
If (itx.Ge.litmax) Exit restart
if (itx >= litmax) exit restart
!
! inner iterations
@@ -410,12 +408,10 @@ Subroutine psb_zrgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist
end if
If (itrace_ > 0) then
if ((mod(itx,itrace_)==0).and.(me == 0))&
& write(*,'(a,i4,3(2x,es10.4))') 'gmres(l): ',itx,errnum,eps*errden
end If
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,eps)
end Do inner
end do inner
if (istop_ == 1) then
x = xt
@@ -435,46 +431,41 @@ Subroutine psb_zrgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist
call psb_geaxpby(zone,w,zone,x,desc_a,info)
end if
End Do restart
If (itrace_ > 0) then
if (me == 0) write(*,'(a,i4,3(2x,es10.4))') 'gmres(l): ',itx,errnum,eps*errden
end If
end do restart
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,1,errnum,errden,eps)
If (Present(err)) then
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).and.(me==0)) Then
write(debug_unit,*) 'rgmres(l) failed to converge to ',eps*errden,&
& ' in ',itx,' iterations '
End If
Deallocate(aux,h,c,s,rs,rst, stat=info)
Call psb_gefree(v,desc_a,info)
Call psb_gefree(w,desc_a,info)
Call psb_gefree(w1,desc_a,info)
Call psb_gefree(xt,desc_a,info)
! restore external global coherence behaviour
call psb_restore_coher(ictxt,isvch)
end if
if (present(iter)) iter = itx
if (errnum > eps*errden) &
& call end_log(methdname,me,itx,errnum,errden,eps)
deallocate(aux,h,c,s,rs,rst, stat=info)
if (info == 0) call psb_gefree(v,desc_a,info)
if (info == 0) call psb_gefree(w,desc_a,info)
if (info == 0) call psb_gefree(w1,desc_a,info)
if (info == 0) call psb_gefree(xt,desc_a,info)
if (info /= 0) then
info=4011
call psb_errpush(info,name)
goto 9999
end if
! restore external global coherence behaviour
call psb_restore_coher(ictxt,isvch)
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
@@ -544,7 +535,7 @@ contains
! .. executable statements ..
!
if( n.le.0 ) return
if( incx.eq.1 .and. incy.eq.1 ) then
if( incx == 1 .and. incy == 1 ) then
!
! code for both increments equal to 1
!