fdurastante/psblas3
1
0
Fork 0
mirror of https://github.com/sfilippone/psblas3.git synced 2026-10-06 22:55:08 +00:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'development' into remap-coarse

Browse Source
This commit is contained in:
sfilippone
2026-06-24 14:02:38 +02:00
parent 66026aca90 d7f3bf536d
commit 4315a2b395
105 changed files with 4728 additions and 921 deletions
Download Patch File Download Diff File Expand all files Collapse all files
base/modules/Makefile
+2 -3
View File
@@ -189,7 +189,7 @@ auxil/psi_acx_mod.o auxil/psi_alcx_mod.o auxil/psi_lcx_mod.o \
auxil/psb_string_mod.o auxil/psb_i2_realloc_mod.o auxil/psb_m_realloc_mod.o auxil/psb_e_realloc_mod.o \
auxil/psb_s_realloc_mod.o \
auxil/psb_d_realloc_mod.o auxil/psb_c_realloc_mod.o auxil/psb_z_realloc_mod.o \
desc/psb_desc_const_mod.o psi_penv_mod.o: psb_const_mod.o
desc/psb_desc_const_mod.o penv/psi_penv_mod.o: psb_const_mod.o
desc/psb_indx_map_mod.o desc/psb_hash_mod.o: psb_realloc_mod.o psb_const_mod.o desc/psb_desc_const_mod.o
@@ -297,8 +297,8 @@ serial/psb_z_csc_mat_mod.o serial/psb_z_csr_mat_mod.o serial/psb_lz_csr_mat_mod.
serial/psb_mat_mod.o: serial/psb_vect_mod.o serial/psb_s_mat_mod.o serial/psb_d_mat_mod.o serial/psb_c_mat_mod.o serial/psb_z_mat_mod.o
serial/psb_serial_mod.o: serial/psb_s_serial_mod.o serial/psb_d_serial_mod.o serial/psb_c_serial_mod.o serial/psb_z_serial_mod.o auxil/psi_serial_mod.o
serial/psb_i2_vect_mod.o: serial/psb_i2_base_vect_mod.o serial/psb_i_base_vect_mod.o
serial/psb_i_vect_mod.o: serial/psb_i_base_vect_mod.o
serial/psb_i2_vect_mod.o: serial/psb_i2_base_vect_mod.o serial/psb_i_vect_mod.o
serial/psb_l_vect_mod.o: serial/psb_l_base_vect_mod.o serial/psb_i_vect_mod.o
serial/psb_s_vect_mod.o: serial/psb_s_base_vect_mod.o serial/psb_i_vect_mod.o
serial/psb_d_vect_mod.o: serial/psb_d_base_vect_mod.o serial/psb_i_vect_mod.o
@@ -447,5 +447,4 @@ clean:
/bin/rm -f $(MODULES) $(OBJS) $(MPFOBJS) *$(.mod)
veryclean: clean
/bin/rm -f *.h
base/modules/desc/psb_desc_mod.F90
+1 -1
View File
@@ -674,7 +674,7 @@ contains
integer(psb_ipk_) :: info
info = 0
if (psb_is_asb_desc(desc)) &
if (desc%is_asb()) &
& call desc%indxmap%set_state(psb_desc_ovl_asb_)
end subroutine psb_cd_set_ovl_asb
base/modules/penv/psi_penv_mod.F90
+8 -8
View File
@@ -1035,6 +1035,8 @@ contains
! !$ call psb_error(ctxt)
! !$ endif
! !$ endif
if (close_) call psb_clear_vect_defaults()
if (close_) call psb_clear_mat_defaults()
#if defined(PSB_SERIAL_MPI)
! Under serial mode, CLOSE has no effect, but reclaim
! the used ctxt number.
@@ -1045,12 +1047,6 @@ contains
else
call psb_close_context(psb_mesg_queue,ctxt)
end if
!if ((ctxt /= mpi_comm_null).and.(ctxt /= mpi_comm_world)) then
if (allocated(ctxt%ctxt)) then
!write(0,*) ctxt%ctxt,mpi_comm_world,mpi_comm_null
if ((ctxt%ctxt /= mpi_comm_world).and.(ctxt%ctxt /= mpi_comm_null)) &
& call mpi_comm_Free(ctxt%ctxt,info)
end if
if (close_) then
if (info == 0) call mpi_op_free(mpi_i2amx_op,info)
if (info == 0) call mpi_op_free(mpi_i2amn_op,info)
@@ -1069,12 +1065,16 @@ contains
if (info == 0) call mpi_op_free(mpi_snrm2_op,info)
if (info == 0) call mpi_op_free(mpi_dnrm2_op,info)
end if
if (allocated(ctxt%ctxt)) then
!!$ write(0,*) ctxt%ctxt,mpi_comm_world,mpi_comm_null,mpi_comm_self
if ((ctxt%ctxt /= mpi_comm_world).and.(ctxt%ctxt /= mpi_comm_null).and.&
& (ctxt%ctxt /= mpi_comm_self)) &
& call mpi_comm_Free(ctxt%ctxt,info)
end if
if (close_) call mpi_finalize(info)
#endif
if (close_) call psb_clear_vect_defaults()
if (close_) call psb_clear_mat_defaults()
end subroutine psb_exit_mpik
base/modules/serial/psb_c_base_vect_mod.F90
+34 -10
View File
@@ -157,6 +157,7 @@ module psb_c_base_vect_mod
procedure, pass(x) :: set_vect => c_base_set_vect
generic, public :: set => set_vect, set_scal
procedure, pass(x) :: get_entry=> c_base_get_entry
procedure, pass(x) :: set_entry=> c_base_set_entry
!
! Gather/scatter. These are needed for MPI interfacing.
! May have to be reworked.
@@ -1275,15 +1276,33 @@ contains
!
function c_base_get_entry(x, index) result(res)
implicit none
class(psb_c_base_vect_type), intent(in) :: x
class(psb_c_base_vect_type), intent(inout) :: x
integer(psb_ipk_), intent(in) :: index
complex(psb_spk_) :: res
res = 0
if (allocated(x%v)) res = x%v(index)
if (allocated(x%v)) then
if (x%is_dev()) call x%sync()
res = x%v(index)
end if
end function c_base_get_entry
subroutine c_base_set_entry(x, index, val)
implicit none
class(psb_c_base_vect_type), intent(inout) :: x
integer(psb_ipk_), intent(in) :: index
complex(psb_spk_) :: val
if (allocated(x%v)) then
if (x%is_dev()) call x%sync()
x%v(index) =val
call x%set_host()
end if
end subroutine c_base_set_entry
!
! Overwrite with absolute value
!
@@ -2677,7 +2696,7 @@ contains
logical, intent(in), optional :: scratch
call psb_realloc(m,n,x%v,info)
call x%asb(m,n,info,scratch)
call x%asb(m,n,info,scratch=scratch)
end subroutine c_base_mlv_bld_n
@@ -2958,23 +2977,26 @@ contains
case(psb_dupl_err_)
do i=1,ncfs
if (any(vv(x%iv(i),:).ne.czero)) then
call psb_errpush(psb_err_duplicate_coo,'vect-asb')
info = psb_err_duplicate_coo
call psb_errpush(info,'mvect-asb')
return
else
vv(x%iv(i),:) = x%v(i,:)
end if
end do
case default
write(psb_err_unit,*) 'Error in vect_asb: unsafe dupl',x%get_dupl()
write(psb_err_unit,*) 'Error in mvect_asb: unsafe dupl',x%get_dupl()
info =-7
end select
call psb_move_alloc(vv,x%v,info)
if (allocated(x%iv)) deallocate(x%iv,stat=info)
else if (x%is_upd().or.x%is_asb().or.scratch_) then
if (x%get_nrows() < m) &
if ((x%get_nrows() < m).or.(x%get_ncols()<n)) &
& call psb_realloc(m,n,x%v,info)
if (info /= 0) &
& call psb_errpush(psb_err_alloc_dealloc_,'vect_asb')
if (info /= 0) then
info = psb_err_alloc_dealloc_
call psb_errpush(psb_err_alloc_dealloc_,'mvect_asb')
end if
else
info = psb_err_invalid_vect_state_
call psb_errpush(info,'vect_asb')
@@ -2982,8 +3004,10 @@ contains
else
if ((x%get_nrows() < m).or.(x%get_ncols()<n)) &
& call psb_realloc(m,n,x%v,info)
if (info /= 0) &
& call psb_errpush(psb_err_alloc_dealloc_,'vect_asb')
if (info /= 0) then
info = psb_err_alloc_dealloc_
call psb_errpush(psb_err_alloc_dealloc_,'mvect_asb')
end if
end if
call x%set_host()
call x%set_asb()
base/modules/serial/psb_c_vect_mod.F90
+16 -5
View File
@@ -108,6 +108,7 @@ module psb_c_vect_mod
procedure, pass(x) :: check_addr => c_vect_check_addr
procedure, pass(x) :: get_entry => c_vect_get_entry
procedure, pass(x) :: set_entry => c_vect_set_entry
procedure, pass(x) :: dot_v => c_vect_dot_v
procedure, pass(x) :: dot_a => c_vect_dot_a
@@ -855,13 +856,21 @@ contains
function c_vect_get_entry(x,index) result(res)
implicit none
class(psb_c_vect_type), intent(in) :: x
class(psb_c_vect_type), intent(inout) :: x
integer(psb_ipk_), intent(in) :: index
complex(psb_spk_) :: res
res = 0
res = czero
if (allocated(x%v)) res = x%v%get_entry(index)
end function c_vect_get_entry
subroutine c_vect_set_entry(x,index,val)
implicit none
class(psb_c_vect_type), intent(inout) :: x
integer(psb_ipk_), intent(in) :: index
complex(psb_spk_) :: val
if (allocated(x%v)) call x%v%set_entry(index,val)
end subroutine c_vect_set_entry
function c_vect_dot_v(n,x,y) result(res)
implicit none
class(psb_c_vect_type), intent(inout) :: x, y
@@ -1660,19 +1669,20 @@ contains
end subroutine c_mvect_bld_x
subroutine c_mvect_bld_n(x,m,n,mold)
subroutine c_mvect_bld_n(x,m,n,mold,scratch)
integer(psb_ipk_), intent(in) :: m,n
class(psb_c_multivect_type), intent(out) :: x
class(psb_c_base_multivect_type), intent(in), optional :: mold
integer(psb_ipk_) :: info
logical, intent(in), optional :: scratch
info = psb_success_
if (present(mold)) then
allocate(x%v,stat=info,mold=mold)
else
allocate(x%v,stat=info, mold=psb_c_get_base_multivect_default())
endif
if (info == psb_success_) call x%v%bld(m,n)
if (info == psb_success_) call x%v%bld(m,n,scratch=scratch)
end subroutine c_mvect_bld_n
@@ -2153,3 +2163,4 @@ contains
!!$ end function c_mvect_asum
end module psb_c_multivect_mod
base/modules/serial/psb_d_base_vect_mod.F90
+38 -10
View File
@@ -157,6 +157,7 @@ module psb_d_base_vect_mod
procedure, pass(x) :: set_vect => d_base_set_vect
generic, public :: set => set_vect, set_scal
procedure, pass(x) :: get_entry=> d_base_get_entry
procedure, pass(x) :: set_entry=> d_base_set_entry
!
! Gather/scatter. These are needed for MPI interfacing.
! May have to be reworked.
@@ -1282,15 +1283,33 @@ contains
!
function d_base_get_entry(x, index) result(res)
implicit none
class(psb_d_base_vect_type), intent(in) :: x
class(psb_d_base_vect_type), intent(inout) :: x
integer(psb_ipk_), intent(in) :: index
real(psb_dpk_) :: res
res = 0
if (allocated(x%v)) res = x%v(index)
if (allocated(x%v)) then
if (x%is_dev()) call x%sync()
res = x%v(index)
end if
end function d_base_get_entry
subroutine d_base_set_entry(x, index, val)
implicit none
class(psb_d_base_vect_type), intent(inout) :: x
integer(psb_ipk_), intent(in) :: index
real(psb_dpk_) :: val
if (allocated(x%v)) then
if (x%is_dev()) call x%sync()
x%v(index) =val
call x%set_host()
end if
end subroutine d_base_set_entry
!
! Overwrite with absolute value
!
@@ -2189,6 +2208,10 @@ contains
res = min(res,abs(x%v(i)))
end do
#else
!
! From M&R&C: if the array is of size zero, MINVAL
! returns the largest positive value
!
res = minval(x%v(1:n))
#endif
end function d_base_min
@@ -2856,7 +2879,7 @@ contains
logical, intent(in), optional :: scratch
call psb_realloc(m,n,x%v,info)
call x%asb(m,n,info,scratch)
call x%asb(m,n,info,scratch=scratch)
end subroutine d_base_mlv_bld_n
@@ -3137,23 +3160,26 @@ contains
case(psb_dupl_err_)
do i=1,ncfs
if (any(vv(x%iv(i),:).ne.dzero)) then
call psb_errpush(psb_err_duplicate_coo,'vect-asb')
info = psb_err_duplicate_coo
call psb_errpush(info,'mvect-asb')
return
else
vv(x%iv(i),:) = x%v(i,:)
end if
end do
case default
write(psb_err_unit,*) 'Error in vect_asb: unsafe dupl',x%get_dupl()
write(psb_err_unit,*) 'Error in mvect_asb: unsafe dupl',x%get_dupl()
info =-7
end select
call psb_move_alloc(vv,x%v,info)
if (allocated(x%iv)) deallocate(x%iv,stat=info)
else if (x%is_upd().or.x%is_asb().or.scratch_) then
if (x%get_nrows() < m) &
if ((x%get_nrows() < m).or.(x%get_ncols()<n)) &
& call psb_realloc(m,n,x%v,info)
if (info /= 0) &
& call psb_errpush(psb_err_alloc_dealloc_,'vect_asb')
if (info /= 0) then
info = psb_err_alloc_dealloc_
call psb_errpush(psb_err_alloc_dealloc_,'mvect_asb')
end if
else
info = psb_err_invalid_vect_state_
call psb_errpush(info,'vect_asb')
@@ -3161,8 +3187,10 @@ contains
else
if ((x%get_nrows() < m).or.(x%get_ncols()<n)) &
& call psb_realloc(m,n,x%v,info)
if (info /= 0) &
& call psb_errpush(psb_err_alloc_dealloc_,'vect_asb')
if (info /= 0) then
info = psb_err_alloc_dealloc_
call psb_errpush(psb_err_alloc_dealloc_,'mvect_asb')
end if
end if
call x%set_host()
call x%set_asb()
base/modules/serial/psb_d_vect_mod.F90
+17 -6
View File
@@ -108,6 +108,7 @@ module psb_d_vect_mod
procedure, pass(x) :: check_addr => d_vect_check_addr
procedure, pass(x) :: get_entry => d_vect_get_entry
procedure, pass(x) :: set_entry => d_vect_set_entry
procedure, pass(x) :: dot_v => d_vect_dot_v
procedure, pass(x) :: dot_a => d_vect_dot_a
@@ -862,13 +863,21 @@ contains
function d_vect_get_entry(x,index) result(res)
implicit none
class(psb_d_vect_type), intent(in) :: x
class(psb_d_vect_type), intent(inout) :: x
integer(psb_ipk_), intent(in) :: index
real(psb_dpk_) :: res
res = 0
res = dzero
if (allocated(x%v)) res = x%v%get_entry(index)
end function d_vect_get_entry
subroutine d_vect_set_entry(x,index,val)
implicit none
class(psb_d_vect_type), intent(inout) :: x
integer(psb_ipk_), intent(in) :: index
real(psb_dpk_) :: val
if (allocated(x%v)) call x%v%set_entry(index,val)
end subroutine d_vect_set_entry
function d_vect_dot_v(n,x,y) result(res)
implicit none
class(psb_d_vect_type), intent(inout) :: x, y
@@ -1430,7 +1439,7 @@ contains
if (allocated(x%v)) then
res = x%v%minreal(n)
else
res = dzero
res = HUGE(done)
end if
end function d_vect_min
@@ -1739,19 +1748,20 @@ contains
end subroutine d_mvect_bld_x
subroutine d_mvect_bld_n(x,m,n,mold)
subroutine d_mvect_bld_n(x,m,n,mold,scratch)
integer(psb_ipk_), intent(in) :: m,n
class(psb_d_multivect_type), intent(out) :: x
class(psb_d_base_multivect_type), intent(in), optional :: mold
integer(psb_ipk_) :: info
logical, intent(in), optional :: scratch
info = psb_success_
if (present(mold)) then
allocate(x%v,stat=info,mold=mold)
else
allocate(x%v,stat=info, mold=psb_d_get_base_multivect_default())
endif
if (info == psb_success_) call x%v%bld(m,n)
if (info == psb_success_) call x%v%bld(m,n,scratch=scratch)
end subroutine d_mvect_bld_n
@@ -2232,3 +2242,4 @@ contains
!!$ end function d_mvect_asum
end module psb_d_multivect_mod
base/modules/serial/psb_i2_base_vect_mod.F90
+13 -8
View File
@@ -1603,7 +1603,7 @@ contains
logical, intent(in), optional :: scratch
call psb_realloc(m,n,x%v,info)
call x%asb(m,n,info,scratch)
call x%asb(m,n,info,scratch=scratch)
end subroutine i2_base_mlv_bld_n
@@ -1884,23 +1884,26 @@ contains
case(psb_dupl_err_)
do i=1,ncfs
if (any(vv(x%iv(i),:).ne.i2zero)) then
call psb_errpush(psb_err_duplicate_coo,'vect-asb')
info = psb_err_duplicate_coo
call psb_errpush(info,'mvect-asb')
return
else
vv(x%iv(i),:) = x%v(i,:)
end if
end do
case default
write(psb_err_unit,*) 'Error in vect_asb: unsafe dupl',x%get_dupl()
write(psb_err_unit,*) 'Error in mvect_asb: unsafe dupl',x%get_dupl()
info =-7
end select
call psb_move_alloc(vv,x%v,info)
if (allocated(x%iv)) deallocate(x%iv,stat=info)
else if (x%is_upd().or.x%is_asb().or.scratch_) then
if (x%get_nrows() < m) &
if ((x%get_nrows() < m).or.(x%get_ncols()<n)) &
& call psb_realloc(m,n,x%v,info)
if (info /= 0) &
& call psb_errpush(psb_err_alloc_dealloc_,'vect_asb')
if (info /= 0) then
info = psb_err_alloc_dealloc_
call psb_errpush(psb_err_alloc_dealloc_,'mvect_asb')
end if
else
info = psb_err_invalid_vect_state_
call psb_errpush(info,'vect_asb')
@@ -1908,8 +1911,10 @@ contains
else
if ((x%get_nrows() < m).or.(x%get_ncols()<n)) &
& call psb_realloc(m,n,x%v,info)
if (info /= 0) &
& call psb_errpush(psb_err_alloc_dealloc_,'vect_asb')
if (info /= 0) then
info = psb_err_alloc_dealloc_
call psb_errpush(psb_err_alloc_dealloc_,'mvect_asb')
end if
end if
call x%set_host()
call x%set_asb()
base/modules/serial/psb_i2_vect_mod.F90
+5 -3
View File
@@ -984,19 +984,20 @@ contains
end subroutine i2_mvect_bld_x
subroutine i2_mvect_bld_n(x,m,n,mold)
subroutine i2_mvect_bld_n(x,m,n,mold,scratch)
integer(psb_ipk_), intent(in) :: m,n
class(psb_i2_multivect_type), intent(out) :: x
class(psb_i2_base_multivect_type), intent(in), optional :: mold
integer(psb_ipk_) :: info
logical, intent(in), optional :: scratch
info = psb_success_
if (present(mold)) then
allocate(x%v,stat=info,mold=mold)
else
allocate(x%v,stat=info, mold=psb_i2_get_base_multivect_default())
endif
if (info == psb_success_) call x%v%bld(m,n)
if (info == psb_success_) call x%v%bld(m,n,scratch=scratch)
end subroutine i2_mvect_bld_n
@@ -1269,3 +1270,4 @@ contains
end module psb_i2_multivect_mod
base/modules/serial/psb_i_base_vect_mod.F90
+13 -8
View File
@@ -1602,7 +1602,7 @@ contains
logical, intent(in), optional :: scratch
call psb_realloc(m,n,x%v,info)
call x%asb(m,n,info,scratch)
call x%asb(m,n,info,scratch=scratch)
end subroutine i_base_mlv_bld_n
@@ -1883,23 +1883,26 @@ contains
case(psb_dupl_err_)
do i=1,ncfs
if (any(vv(x%iv(i),:).ne.izero)) then
call psb_errpush(psb_err_duplicate_coo,'vect-asb')
info = psb_err_duplicate_coo
call psb_errpush(info,'mvect-asb')
return
else
vv(x%iv(i),:) = x%v(i,:)
end if
end do
case default
write(psb_err_unit,*) 'Error in vect_asb: unsafe dupl',x%get_dupl()
write(psb_err_unit,*) 'Error in mvect_asb: unsafe dupl',x%get_dupl()
info =-7
end select
call psb_move_alloc(vv,x%v,info)
if (allocated(x%iv)) deallocate(x%iv,stat=info)
else if (x%is_upd().or.x%is_asb().or.scratch_) then
if (x%get_nrows() < m) &
if ((x%get_nrows() < m).or.(x%get_ncols()<n)) &
& call psb_realloc(m,n,x%v,info)
if (info /= 0) &
& call psb_errpush(psb_err_alloc_dealloc_,'vect_asb')
if (info /= 0) then
info = psb_err_alloc_dealloc_
call psb_errpush(psb_err_alloc_dealloc_,'mvect_asb')
end if
else
info = psb_err_invalid_vect_state_
call psb_errpush(info,'vect_asb')
@@ -1907,8 +1910,10 @@ contains
else
if ((x%get_nrows() < m).or.(x%get_ncols()<n)) &
& call psb_realloc(m,n,x%v,info)
if (info /= 0) &
& call psb_errpush(psb_err_alloc_dealloc_,'vect_asb')
if (info /= 0) then
info = psb_err_alloc_dealloc_
call psb_errpush(psb_err_alloc_dealloc_,'mvect_asb')
end if
end if
call x%set_host()
call x%set_asb()
base/modules/serial/psb_i_vect_mod.F90
+5 -3
View File
@@ -983,19 +983,20 @@ contains
end subroutine i_mvect_bld_x
subroutine i_mvect_bld_n(x,m,n,mold)
subroutine i_mvect_bld_n(x,m,n,mold,scratch)
integer(psb_ipk_), intent(in) :: m,n
class(psb_i_multivect_type), intent(out) :: x
class(psb_i_base_multivect_type), intent(in), optional :: mold
integer(psb_ipk_) :: info
logical, intent(in), optional :: scratch
info = psb_success_
if (present(mold)) then
allocate(x%v,stat=info,mold=mold)
else
allocate(x%v,stat=info, mold=psb_i_get_base_multivect_default())
endif
if (info == psb_success_) call x%v%bld(m,n)
if (info == psb_success_) call x%v%bld(m,n,scratch=scratch)
end subroutine i_mvect_bld_n
@@ -1268,3 +1269,4 @@ contains
end module psb_i_multivect_mod
base/modules/serial/psb_l_base_vect_mod.F90
+13 -8
View File
@@ -1603,7 +1603,7 @@ contains
logical, intent(in), optional :: scratch
call psb_realloc(m,n,x%v,info)
call x%asb(m,n,info,scratch)
call x%asb(m,n,info,scratch=scratch)
end subroutine l_base_mlv_bld_n
@@ -1884,23 +1884,26 @@ contains
case(psb_dupl_err_)
do i=1,ncfs
if (any(vv(x%iv(i),:).ne.lzero)) then
call psb_errpush(psb_err_duplicate_coo,'vect-asb')
info = psb_err_duplicate_coo
call psb_errpush(info,'mvect-asb')
return
else
vv(x%iv(i),:) = x%v(i,:)
end if
end do
case default
write(psb_err_unit,*) 'Error in vect_asb: unsafe dupl',x%get_dupl()
write(psb_err_unit,*) 'Error in mvect_asb: unsafe dupl',x%get_dupl()
info =-7
end select
call psb_move_alloc(vv,x%v,info)
if (allocated(x%iv)) deallocate(x%iv,stat=info)
else if (x%is_upd().or.x%is_asb().or.scratch_) then
if (x%get_nrows() < m) &
if ((x%get_nrows() < m).or.(x%get_ncols()<n)) &
& call psb_realloc(m,n,x%v,info)
if (info /= 0) &
& call psb_errpush(psb_err_alloc_dealloc_,'vect_asb')
if (info /= 0) then
info = psb_err_alloc_dealloc_
call psb_errpush(psb_err_alloc_dealloc_,'mvect_asb')
end if
else
info = psb_err_invalid_vect_state_
call psb_errpush(info,'vect_asb')
@@ -1908,8 +1911,10 @@ contains
else
if ((x%get_nrows() < m).or.(x%get_ncols()<n)) &
& call psb_realloc(m,n,x%v,info)
if (info /= 0) &
& call psb_errpush(psb_err_alloc_dealloc_,'vect_asb')
if (info /= 0) then
info = psb_err_alloc_dealloc_
call psb_errpush(psb_err_alloc_dealloc_,'mvect_asb')
end if
end if
call x%set_host()
call x%set_asb()
base/modules/serial/psb_l_vect_mod.F90
+5 -3
View File
@@ -984,19 +984,20 @@ contains
end subroutine l_mvect_bld_x
subroutine l_mvect_bld_n(x,m,n,mold)
subroutine l_mvect_bld_n(x,m,n,mold,scratch)
integer(psb_ipk_), intent(in) :: m,n
class(psb_l_multivect_type), intent(out) :: x
class(psb_l_base_multivect_type), intent(in), optional :: mold
integer(psb_ipk_) :: info
logical, intent(in), optional :: scratch
info = psb_success_
if (present(mold)) then
allocate(x%v,stat=info,mold=mold)
else
allocate(x%v,stat=info, mold=psb_l_get_base_multivect_default())
endif
if (info == psb_success_) call x%v%bld(m,n)
if (info == psb_success_) call x%v%bld(m,n,scratch=scratch)
end subroutine l_mvect_bld_n
@@ -1269,3 +1270,4 @@ contains
end module psb_l_multivect_mod
base/modules/serial/psb_s_base_vect_mod.F90
+38 -10
View File
@@ -157,6 +157,7 @@ module psb_s_base_vect_mod
procedure, pass(x) :: set_vect => s_base_set_vect
generic, public :: set => set_vect, set_scal
procedure, pass(x) :: get_entry=> s_base_get_entry
procedure, pass(x) :: set_entry=> s_base_set_entry
!
! Gather/scatter. These are needed for MPI interfacing.
! May have to be reworked.
@@ -1282,15 +1283,33 @@ contains
!
function s_base_get_entry(x, index) result(res)
implicit none
class(psb_s_base_vect_type), intent(in) :: x
class(psb_s_base_vect_type), intent(inout) :: x
integer(psb_ipk_), intent(in) :: index
real(psb_spk_) :: res
res = 0
if (allocated(x%v)) res = x%v(index)
if (allocated(x%v)) then
if (x%is_dev()) call x%sync()
res = x%v(index)
end if
end function s_base_get_entry
subroutine s_base_set_entry(x, index, val)
implicit none
class(psb_s_base_vect_type), intent(inout) :: x
integer(psb_ipk_), intent(in) :: index
real(psb_spk_) :: val
if (allocated(x%v)) then
if (x%is_dev()) call x%sync()
x%v(index) =val
call x%set_host()
end if
end subroutine s_base_set_entry
!
! Overwrite with absolute value
!
@@ -2189,6 +2208,10 @@ contains
res = min(res,abs(x%v(i)))
end do
#else
!
! From M&R&C: if the array is of size zero, MINVAL
! returns the largest positive value
!
res = minval(x%v(1:n))
#endif
end function s_base_min
@@ -2856,7 +2879,7 @@ contains
logical, intent(in), optional :: scratch
call psb_realloc(m,n,x%v,info)
call x%asb(m,n,info,scratch)
call x%asb(m,n,info,scratch=scratch)
end subroutine s_base_mlv_bld_n
@@ -3137,23 +3160,26 @@ contains
case(psb_dupl_err_)
do i=1,ncfs
if (any(vv(x%iv(i),:).ne.szero)) then
call psb_errpush(psb_err_duplicate_coo,'vect-asb')
info = psb_err_duplicate_coo
call psb_errpush(info,'mvect-asb')
return
else
vv(x%iv(i),:) = x%v(i,:)
end if
end do
case default
write(psb_err_unit,*) 'Error in vect_asb: unsafe dupl',x%get_dupl()
write(psb_err_unit,*) 'Error in mvect_asb: unsafe dupl',x%get_dupl()
info =-7
end select
call psb_move_alloc(vv,x%v,info)
if (allocated(x%iv)) deallocate(x%iv,stat=info)
else if (x%is_upd().or.x%is_asb().or.scratch_) then
if (x%get_nrows() < m) &
if ((x%get_nrows() < m).or.(x%get_ncols()<n)) &
& call psb_realloc(m,n,x%v,info)
if (info /= 0) &
& call psb_errpush(psb_err_alloc_dealloc_,'vect_asb')
if (info /= 0) then
info = psb_err_alloc_dealloc_
call psb_errpush(psb_err_alloc_dealloc_,'mvect_asb')
end if
else
info = psb_err_invalid_vect_state_
call psb_errpush(info,'vect_asb')
@@ -3161,8 +3187,10 @@ contains
else
if ((x%get_nrows() < m).or.(x%get_ncols()<n)) &
& call psb_realloc(m,n,x%v,info)
if (info /= 0) &
& call psb_errpush(psb_err_alloc_dealloc_,'vect_asb')
if (info /= 0) then
info = psb_err_alloc_dealloc_
call psb_errpush(psb_err_alloc_dealloc_,'mvect_asb')
end if
end if
call x%set_host()
call x%set_asb()
base/modules/serial/psb_s_vect_mod.F90
+17 -6
View File
@@ -108,6 +108,7 @@ module psb_s_vect_mod
procedure, pass(x) :: check_addr => s_vect_check_addr
procedure, pass(x) :: get_entry => s_vect_get_entry
procedure, pass(x) :: set_entry => s_vect_set_entry
procedure, pass(x) :: dot_v => s_vect_dot_v
procedure, pass(x) :: dot_a => s_vect_dot_a
@@ -862,13 +863,21 @@ contains
function s_vect_get_entry(x,index) result(res)
implicit none
class(psb_s_vect_type), intent(in) :: x
class(psb_s_vect_type), intent(inout) :: x
integer(psb_ipk_), intent(in) :: index
real(psb_spk_) :: res
res = 0
res = szero
if (allocated(x%v)) res = x%v%get_entry(index)
end function s_vect_get_entry
subroutine s_vect_set_entry(x,index,val)
implicit none
class(psb_s_vect_type), intent(inout) :: x
integer(psb_ipk_), intent(in) :: index
real(psb_spk_) :: val
if (allocated(x%v)) call x%v%set_entry(index,val)
end subroutine s_vect_set_entry
function s_vect_dot_v(n,x,y) result(res)
implicit none
class(psb_s_vect_type), intent(inout) :: x, y
@@ -1430,7 +1439,7 @@ contains
if (allocated(x%v)) then
res = x%v%minreal(n)
else
res = szero
res = HUGE(sone)
end if
end function s_vect_min
@@ -1739,19 +1748,20 @@ contains
end subroutine s_mvect_bld_x
subroutine s_mvect_bld_n(x,m,n,mold)
subroutine s_mvect_bld_n(x,m,n,mold,scratch)
integer(psb_ipk_), intent(in) :: m,n
class(psb_s_multivect_type), intent(out) :: x
class(psb_s_base_multivect_type), intent(in), optional :: mold
integer(psb_ipk_) :: info
logical, intent(in), optional :: scratch
info = psb_success_
if (present(mold)) then
allocate(x%v,stat=info,mold=mold)
else
allocate(x%v,stat=info, mold=psb_s_get_base_multivect_default())
endif
if (info == psb_success_) call x%v%bld(m,n)
if (info == psb_success_) call x%v%bld(m,n,scratch=scratch)
end subroutine s_mvect_bld_n
@@ -2232,3 +2242,4 @@ contains
!!$ end function s_mvect_asum
end module psb_s_multivect_mod
base/modules/serial/psb_z_base_vect_mod.F90
+34 -10
View File
@@ -157,6 +157,7 @@ module psb_z_base_vect_mod
procedure, pass(x) :: set_vect => z_base_set_vect
generic, public :: set => set_vect, set_scal
procedure, pass(x) :: get_entry=> z_base_get_entry
procedure, pass(x) :: set_entry=> z_base_set_entry
!
! Gather/scatter. These are needed for MPI interfacing.
! May have to be reworked.
@@ -1275,15 +1276,33 @@ contains
!
function z_base_get_entry(x, index) result(res)
implicit none
class(psb_z_base_vect_type), intent(in) :: x
class(psb_z_base_vect_type), intent(inout) :: x
integer(psb_ipk_), intent(in) :: index
complex(psb_dpk_) :: res
res = 0
if (allocated(x%v)) res = x%v(index)
if (allocated(x%v)) then
if (x%is_dev()) call x%sync()
res = x%v(index)
end if
end function z_base_get_entry
subroutine z_base_set_entry(x, index, val)
implicit none
class(psb_z_base_vect_type), intent(inout) :: x
integer(psb_ipk_), intent(in) :: index
complex(psb_dpk_) :: val
if (allocated(x%v)) then
if (x%is_dev()) call x%sync()
x%v(index) =val
call x%set_host()
end if
end subroutine z_base_set_entry
!
! Overwrite with absolute value
!
@@ -2677,7 +2696,7 @@ contains
logical, intent(in), optional :: scratch
call psb_realloc(m,n,x%v,info)
call x%asb(m,n,info,scratch)
call x%asb(m,n,info,scratch=scratch)
end subroutine z_base_mlv_bld_n
@@ -2958,23 +2977,26 @@ contains
case(psb_dupl_err_)
do i=1,ncfs
if (any(vv(x%iv(i),:).ne.zzero)) then
call psb_errpush(psb_err_duplicate_coo,'vect-asb')
info = psb_err_duplicate_coo
call psb_errpush(info,'mvect-asb')
return
else
vv(x%iv(i),:) = x%v(i,:)
end if
end do
case default
write(psb_err_unit,*) 'Error in vect_asb: unsafe dupl',x%get_dupl()
write(psb_err_unit,*) 'Error in mvect_asb: unsafe dupl',x%get_dupl()
info =-7
end select
call psb_move_alloc(vv,x%v,info)
if (allocated(x%iv)) deallocate(x%iv,stat=info)
else if (x%is_upd().or.x%is_asb().or.scratch_) then
if (x%get_nrows() < m) &
if ((x%get_nrows() < m).or.(x%get_ncols()<n)) &
& call psb_realloc(m,n,x%v,info)
if (info /= 0) &
& call psb_errpush(psb_err_alloc_dealloc_,'vect_asb')
if (info /= 0) then
info = psb_err_alloc_dealloc_
call psb_errpush(psb_err_alloc_dealloc_,'mvect_asb')
end if
else
info = psb_err_invalid_vect_state_
call psb_errpush(info,'vect_asb')
@@ -2982,8 +3004,10 @@ contains
else
if ((x%get_nrows() < m).or.(x%get_ncols()<n)) &
& call psb_realloc(m,n,x%v,info)
if (info /= 0) &
& call psb_errpush(psb_err_alloc_dealloc_,'vect_asb')
if (info /= 0) then
info = psb_err_alloc_dealloc_
call psb_errpush(psb_err_alloc_dealloc_,'mvect_asb')
end if
end if
call x%set_host()
call x%set_asb()
base/modules/serial/psb_z_vect_mod.F90
+16 -5
View File
@@ -108,6 +108,7 @@ module psb_z_vect_mod
procedure, pass(x) :: check_addr => z_vect_check_addr
procedure, pass(x) :: get_entry => z_vect_get_entry
procedure, pass(x) :: set_entry => z_vect_set_entry
procedure, pass(x) :: dot_v => z_vect_dot_v
procedure, pass(x) :: dot_a => z_vect_dot_a
@@ -855,13 +856,21 @@ contains
function z_vect_get_entry(x,index) result(res)
implicit none
class(psb_z_vect_type), intent(in) :: x
class(psb_z_vect_type), intent(inout) :: x
integer(psb_ipk_), intent(in) :: index
complex(psb_dpk_) :: res
res = 0
res = zzero
if (allocated(x%v)) res = x%v%get_entry(index)
end function z_vect_get_entry
subroutine z_vect_set_entry(x,index,val)
implicit none
class(psb_z_vect_type), intent(inout) :: x
integer(psb_ipk_), intent(in) :: index
complex(psb_dpk_) :: val
if (allocated(x%v)) call x%v%set_entry(index,val)
end subroutine z_vect_set_entry
function z_vect_dot_v(n,x,y) result(res)
implicit none
class(psb_z_vect_type), intent(inout) :: x, y
@@ -1660,19 +1669,20 @@ contains
end subroutine z_mvect_bld_x
subroutine z_mvect_bld_n(x,m,n,mold)
subroutine z_mvect_bld_n(x,m,n,mold,scratch)
integer(psb_ipk_), intent(in) :: m,n
class(psb_z_multivect_type), intent(out) :: x
class(psb_z_base_multivect_type), intent(in), optional :: mold
integer(psb_ipk_) :: info
logical, intent(in), optional :: scratch
info = psb_success_
if (present(mold)) then
allocate(x%v,stat=info,mold=mold)
else
allocate(x%v,stat=info, mold=psb_z_get_base_multivect_default())
endif
if (info == psb_success_) call x%v%bld(m,n)
if (info == psb_success_) call x%v%bld(m,n,scratch=scratch)
end subroutine z_mvect_bld_n
@@ -2153,3 +2163,4 @@ contains
!!$ end function z_mvect_asum
end module psb_z_multivect_mod
base/psblas/psb_damax.f90
+2 -2
View File
@@ -689,11 +689,11 @@ function psb_dmin_vect(x, desc_a, info,global) result(res)
goto 9999
end if
! compute local max
! compute local min
if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then
res = x%minreal(desc_a%get_local_rows())
else
res = dzero
res = HUGE(done)
end if
! compute global min
base/psblas/psb_samax.f90
+2 -2
View File
@@ -689,11 +689,11 @@ function psb_smin_vect(x, desc_a, info,global) result(res)
goto 9999
end if
! compute local max
! compute local min
if ((desc_a%get_local_rows() > 0).and.(m /= 0)) then
res = x%minreal(desc_a%get_local_rows())
else
res = szero
res = HUGE(sone)
end if
! compute global min
base/tools/psb_casb.f90
+1 -1
View File
@@ -309,7 +309,7 @@ subroutine psb_casb_multivect(x, desc_a, info, mold, scratch,n)
dupl_ = x%get_dupl()
if (scratch_) then
call x%free(info)
call x%bld(ncol,n_,mold=mold)
call x%bld(ncol,n_,mold=mold,scratch=.true.)
else
call x%asb(ncol,n_,info)
! ..update halo elements..
base/tools/psb_dasb.f90
+1 -1
View File
@@ -309,7 +309,7 @@ subroutine psb_dasb_multivect(x, desc_a, info, mold, scratch,n)
dupl_ = x%get_dupl()
if (scratch_) then
call x%free(info)
call x%bld(ncol,n_,mold=mold)
call x%bld(ncol,n_,mold=mold,scratch=.true.)
else
call x%asb(ncol,n_,info)
! ..update halo elements..
base/tools/psb_iasb.f90
+1 -1
View File
@@ -309,7 +309,7 @@ subroutine psb_iasb_multivect(x, desc_a, info, mold, scratch,n)
dupl_ = x%get_dupl()
if (scratch_) then
call x%free(info)
call x%bld(ncol,n_,mold=mold)
call x%bld(ncol,n_,mold=mold,scratch=.true.)
else
call x%asb(ncol,n_,info)
! ..update halo elements..
base/tools/psb_lasb.f90
+1 -1
View File
@@ -309,7 +309,7 @@ subroutine psb_lasb_multivect(x, desc_a, info, mold, scratch,n)
dupl_ = x%get_dupl()
if (scratch_) then
call x%free(info)
call x%bld(ncol,n_,mold=mold)
call x%bld(ncol,n_,mold=mold,scratch=.true.)
else
call x%asb(ncol,n_,info)
! ..update halo elements..
base/tools/psb_sasb.f90
+1 -1
View File
@@ -309,7 +309,7 @@ subroutine psb_sasb_multivect(x, desc_a, info, mold, scratch,n)
dupl_ = x%get_dupl()
if (scratch_) then
call x%free(info)
call x%bld(ncol,n_,mold=mold)
call x%bld(ncol,n_,mold=mold,scratch=.true.)
else
call x%asb(ncol,n_,info)
! ..update halo elements..
base/tools/psb_zasb.f90
+1 -1
View File
@@ -309,7 +309,7 @@ subroutine psb_zasb_multivect(x, desc_a, info, mold, scratch,n)
dupl_ = x%get_dupl()
if (scratch_) then
call x%free(info)
call x%bld(ncol,n_,mold=mold)
call x%bld(ncol,n_,mold=mold,scratch=.true.)
else
call x%asb(ncol,n_,info)
! ..update halo elements..
cbind/base/psb_base_cbind_mod.f90
+4
View File
@@ -1,6 +1,10 @@
module psb_base_cbind_mod
use psb_objhandle_mod
use psb_cpenv_mod
use psb_s_serial_cbind_mod
use psb_d_serial_cbind_mod
use psb_c_serial_cbind_mod
use psb_z_serial_cbind_mod
use psb_base_tools_cbind_mod
use psb_s_tools_cbind_mod
use psb_d_tools_cbind_mod
cbind/base/psb_c_cbase.h
+4 -2
View File
@@ -64,8 +64,10 @@ psb_i_t psb_c_cspasb_opt(psb_c_cspmat *mh, psb_c_descriptor *cdh,
const char *afmt, psb_i_t upd, psb_i_t dupl);
psb_i_t psb_c_csprn(psb_c_cspmat *mh, psb_c_descriptor *cdh, _Bool clear);
psb_i_t psb_c_cmat_name_print(psb_c_cspmat *mh, char *name);
psb_i_t psb_c_cvect_set_scal(psb_c_cvector *xh, psb_c_t val);
psb_i_t psb_c_cvect_set_vect(psb_c_cvector *xh, psb_c_t *val, psb_i_t n);
psb_i_t psb_c_cvect_set_scal(psb_c_cvector *xh, psb_c_t val);
psb_i_t psb_c_cvect_set_vect(psb_c_cvector *xh, psb_c_t *val, psb_i_t n);
psb_c_t psb_c_cvect_get_entry(psb_c_cvector *xh, psb_i_t index);
psb_i_t psb_c_cvect_set_entry(psb_c_cvector *xh, psb_i_t index, psb_c_t val);
/* psblas computational routines */
psb_c_t psb_c_cgedot(psb_c_cvector *xh, psb_c_cvector *yh, psb_c_descriptor *cdh);
cbind/base/psb_c_dbase.h
+4 -2
View File
@@ -65,8 +65,10 @@ psb_i_t psb_c_dspasb_opt(psb_c_dspmat *mh, psb_c_descriptor *cdh,
const char *afmt, psb_i_t upd, psb_i_t dupl);
psb_i_t psb_c_dsprn(psb_c_dspmat *mh, psb_c_descriptor *cdh, _Bool clear);
psb_i_t psb_c_dmat_name_print(psb_c_dspmat *mh, char *name);
psb_i_t psb_c_dvect_set_scal(psb_c_dvector *xh, psb_d_t val);
psb_i_t psb_c_dvect_set_vect(psb_c_dvector *xh, psb_d_t *val, psb_i_t n);
psb_i_t psb_c_dvect_set_scal(psb_c_dvector *xh, psb_d_t val);
psb_i_t psb_c_dvect_set_vect(psb_c_dvector *xh, psb_d_t *val, psb_i_t n);
psb_d_t psb_c_dvect_get_entry(psb_c_dvector *xh, psb_i_t index);
psb_i_t psb_c_dvect_set_entry(psb_c_dvector *xh, psb_i_t index, psb_d_t val);
/* psblas computational routines */
psb_d_t psb_c_dgedot(psb_c_dvector *xh, psb_c_dvector *yh, psb_c_descriptor *cdh);
cbind/base/psb_c_psblas_cbind_mod.f90
+1 -2
View File
@@ -429,14 +429,13 @@ contains
function psb_c_cgecmp(xh,ch,zh,cdh) bind(c) result(res)
implicit none
integer(psb_c_ipk_) :: res
real(c_float_complex), value :: ch
type(psb_c_cvector) :: xh,zh
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_c_vect_type), pointer :: xp,zp
integer(psb_c_ipk_) :: info
real(c_float_complex), value :: ch
res = -1
cbind/base/psb_c_sbase.h
+5 -2
View File
@@ -65,8 +65,11 @@ psb_i_t psb_c_scopy_mat(psb_c_sspmat *ah,psb_c_sspmat *bh,psb_c_descriptor *cd
const char *afmt, psb_i_t upd, psb_i_t dupl);
psb_i_t psb_c_ssprn(psb_c_sspmat *mh, psb_c_descriptor *cdh, _Bool clear);
psb_i_t psb_c_smat_name_print(psb_c_sspmat *mh, char *name);
psb_i_t psb_c_svect_set_scal(psb_c_svector *xh, psb_s_t val);
psb_i_t psb_c_svect_set_vect(psb_c_svector *xh, psb_s_t *val, psb_i_t n);
psb_i_t psb_c_svect_set_scal(psb_c_svector *xh, psb_s_t val);
psb_i_t psb_c_svect_set_vect(psb_c_svector *xh, psb_s_t *val, psb_i_t n);
psb_s_t psb_c_svect_get_entry(psb_c_svector *xh, psb_i_t index);
psb_i_t psb_c_svect_set_entry(psb_c_svector *xh, psb_i_t index, psb_s_t val);
/* psblas computational routines */
psb_s_t psb_c_sgedot(psb_c_svector *xh, psb_c_svector *yh, psb_c_descriptor *cdh);
cbind/base/psb_c_serial_cbind_mod.F90
+45
View File
@@ -200,6 +200,51 @@ contains
end function psb_c_cvect_set_vect
function psb_c_cvect_set_entry(x,index,val) bind(c) result(info)
use psb_base_mod
implicit none
type(psb_c_cvector) :: x
type(psb_c_vect_type), pointer :: xp
integer(psb_c_ipk_) :: info
integer(psb_c_ipk_), value :: index
complex(c_float_complex), value :: val
integer(psb_c_ipk_) :: ixb
info = -1;
if (c_associated(x%item)) then
call c_f_pointer(x%item,xp)
else
return
end if
ixb = psb_c_get_index_base()
call xp%set_entry((index+(1-ixb)),val)
info = 0
end function psb_c_cvect_set_entry
function psb_c_cvect_get_entry(x,index) bind(c) result(res)
use psb_base_mod
implicit none
type(psb_c_cvector) :: x
type(psb_c_vect_type), pointer :: xp
integer(psb_c_ipk_), value :: index
complex(c_float_complex) :: res
integer(psb_c_ipk_) :: ixb
if (c_associated(x%item)) then
call c_f_pointer(x%item,xp)
else
return
end if
ixb = psb_c_get_index_base()
res = xp%get_entry((index+(1-ixb)))
end function psb_c_cvect_get_entry
function psb_c_cvect_clone(xh,yh) bind(c) result(info)
implicit none
cbind/base/psb_c_tools_cbind_mod.F90
+11 -5
View File
@@ -281,13 +281,19 @@ contains
end if
ixb = psb_c_get_index_base()
if (ixb == 1) then
call psb_geins(nz,irw(1:nz),val(1:nz),&
& xp,descp,info)
else
select case(ixb)
case (0)
!write(0,*) 'C_GEINS: IDX_BASE',ixb,' :',irw(1:nz)+(1-ixb),val(1:nz)
call psb_geins(nz,(irw(1:nz)+(1-ixb)),val(1:nz),&
& xp,descp,info)
end if
case(1)
!write(0,*) 'C_GEINS: IDX_BASE',ixb,' :',irw(1:nz),val(1:nz)
call psb_geins(nz,irw(1:nz),val(1:nz),&
& xp,descp,info)
case default
write(0,*) 'C_GEINS: Unkonwn inndex base ',ixb
info =-2
end select
res = min(0,info)
cbind/base/psb_c_zbase.h
+4 -2
View File
@@ -66,8 +66,10 @@ psb_i_t psb_c_zspasb_opt(psb_c_zspmat *mh, psb_c_descriptor *cdh,
const char *afmt, psb_i_t upd, psb_i_t dupl);
psb_i_t psb_c_zsprn(psb_c_zspmat *mh, psb_c_descriptor *cdh, _Bool clear);
psb_i_t psb_c_zmat_name_print(psb_c_zspmat *mh, char *name);
psb_i_t psb_c_zvect_set_scal(psb_c_zvector *xh, psb_z_t val);
psb_i_t psb_c_zvect_set_vect(psb_c_zvector *xh, psb_z_t *val, psb_i_t n);
psb_i_t psb_c_zvect_set_scal(psb_c_zvector *xh, psb_z_t val);
psb_i_t psb_c_zvect_set_vect(psb_c_zvector *xh, psb_z_t *val, psb_i_t n);
psb_z_t psb_c_zvect_get_entry(psb_c_zvector *xh, psb_i_t index);
psb_i_t psb_c_zvect_set_entry(psb_c_zvector *xh, psb_i_t index, psb_z_t val);
/* psblas computational routines */
psb_z_t psb_c_zgedot(psb_c_zvector *xh, psb_c_zvector *yh, psb_c_descriptor *cdh);
cbind/base/psb_d_psblas_cbind_mod.f90
+1 -2
View File
@@ -429,14 +429,13 @@ contains
function psb_c_dgecmp(xh,ch,zh,cdh) bind(c) result(res)
implicit none
integer(psb_c_ipk_) :: res
real(c_double), value :: ch
type(psb_c_dvector) :: xh,zh
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_d_vect_type), pointer :: xp,zp
integer(psb_c_ipk_) :: info
real(c_double), value :: ch
res = -1
cbind/base/psb_d_serial_cbind_mod.F90
+45
View File
@@ -200,6 +200,51 @@ contains
end function psb_c_dvect_set_vect
function psb_c_dvect_set_entry(x,index,val) bind(c) result(info)
use psb_base_mod
implicit none
type(psb_c_dvector) :: x
type(psb_d_vect_type), pointer :: xp
integer(psb_c_ipk_) :: info
integer(psb_c_ipk_), value :: index
real(c_double), value :: val
integer(psb_c_ipk_) :: ixb
info = -1;
if (c_associated(x%item)) then
call c_f_pointer(x%item,xp)
else
return
end if
ixb = psb_c_get_index_base()
call xp%set_entry((index+(1-ixb)),val)
info = 0
end function psb_c_dvect_set_entry
function psb_c_dvect_get_entry(x,index) bind(c) result(res)
use psb_base_mod
implicit none
type(psb_c_dvector) :: x
type(psb_d_vect_type), pointer :: xp
integer(psb_c_ipk_), value :: index
real(c_double) :: res
integer(psb_c_ipk_) :: ixb
if (c_associated(x%item)) then
call c_f_pointer(x%item,xp)
else
return
end if
ixb = psb_c_get_index_base()
res = xp%get_entry((index+(1-ixb)))
end function psb_c_dvect_get_entry
function psb_c_dvect_clone(xh,yh) bind(c) result(info)
implicit none
cbind/base/psb_d_tools_cbind_mod.F90
+11 -5
View File
@@ -281,13 +281,19 @@ contains
end if
ixb = psb_c_get_index_base()
if (ixb == 1) then
call psb_geins(nz,irw(1:nz),val(1:nz),&
& xp,descp,info)
else
select case(ixb)
case (0)
!write(0,*) 'C_GEINS: IDX_BASE',ixb,' :',irw(1:nz)+(1-ixb),val(1:nz)
call psb_geins(nz,(irw(1:nz)+(1-ixb)),val(1:nz),&
& xp,descp,info)
end if
case(1)
!write(0,*) 'C_GEINS: IDX_BASE',ixb,' :',irw(1:nz),val(1:nz)
call psb_geins(nz,irw(1:nz),val(1:nz),&
& xp,descp,info)
case default
write(0,*) 'C_GEINS: Unkonwn inndex base ',ixb
info =-2
end select
res = min(0,info)
cbind/base/psb_s_psblas_cbind_mod.f90
+1 -2
View File
@@ -429,14 +429,13 @@ contains
function psb_c_sgecmp(xh,ch,zh,cdh) bind(c) result(res)
implicit none
integer(psb_c_ipk_) :: res
real(c_float), value :: ch
type(psb_c_svector) :: xh,zh
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_s_vect_type), pointer :: xp,zp
integer(psb_c_ipk_) :: info
real(c_float), value :: ch
res = -1
cbind/base/psb_s_serial_cbind_mod.F90
+45
View File
@@ -200,6 +200,51 @@ contains
end function psb_c_svect_set_vect
function psb_c_svect_set_entry(x,index,val) bind(c) result(info)
use psb_base_mod
implicit none
type(psb_c_svector) :: x
type(psb_s_vect_type), pointer :: xp
integer(psb_c_ipk_) :: info
integer(psb_c_ipk_), value :: index
real(c_float), value :: val
integer(psb_c_ipk_) :: ixb
info = -1;
if (c_associated(x%item)) then
call c_f_pointer(x%item,xp)
else
return
end if
ixb = psb_c_get_index_base()
call xp%set_entry((index+(1-ixb)),val)
info = 0
end function psb_c_svect_set_entry
function psb_c_svect_get_entry(x,index) bind(c) result(res)
use psb_base_mod
implicit none
type(psb_c_svector) :: x
type(psb_s_vect_type), pointer :: xp
integer(psb_c_ipk_), value :: index
real(c_float) :: res
integer(psb_c_ipk_) :: ixb
if (c_associated(x%item)) then
call c_f_pointer(x%item,xp)
else
return
end if
ixb = psb_c_get_index_base()
res = xp%get_entry((index+(1-ixb)))
end function psb_c_svect_get_entry
function psb_c_svect_clone(xh,yh) bind(c) result(info)
implicit none
cbind/base/psb_s_tools_cbind_mod.F90
+11 -5
View File
@@ -281,13 +281,19 @@ contains
end if
ixb = psb_c_get_index_base()
if (ixb == 1) then
call psb_geins(nz,irw(1:nz),val(1:nz),&
& xp,descp,info)
else
select case(ixb)
case (0)
!write(0,*) 'C_GEINS: IDX_BASE',ixb,' :',irw(1:nz)+(1-ixb),val(1:nz)
call psb_geins(nz,(irw(1:nz)+(1-ixb)),val(1:nz),&
& xp,descp,info)
end if
case(1)
!write(0,*) 'C_GEINS: IDX_BASE',ixb,' :',irw(1:nz),val(1:nz)
call psb_geins(nz,irw(1:nz),val(1:nz),&
& xp,descp,info)
case default
write(0,*) 'C_GEINS: Unkonwn inndex base ',ixb
info =-2
end select
res = min(0,info)
cbind/base/psb_z_psblas_cbind_mod.f90
+1 -2
View File
@@ -429,14 +429,13 @@ contains
function psb_c_zgecmp(xh,ch,zh,cdh) bind(c) result(res)
implicit none
integer(psb_c_ipk_) :: res
real(c_double_complex), value :: ch
type(psb_c_zvector) :: xh,zh
type(psb_c_descriptor) :: cdh
type(psb_desc_type), pointer :: descp
type(psb_z_vect_type), pointer :: xp,zp
integer(psb_c_ipk_) :: info
real(c_double_complex), value :: ch
res = -1
cbind/base/psb_z_serial_cbind_mod.F90
+45
View File
@@ -200,6 +200,51 @@ contains
end function psb_c_zvect_set_vect
function psb_c_zvect_set_entry(x,index,val) bind(c) result(info)
use psb_base_mod
implicit none
type(psb_c_zvector) :: x
type(psb_z_vect_type), pointer :: xp
integer(psb_c_ipk_) :: info
integer(psb_c_ipk_), value :: index
complex(c_double_complex), value :: val
integer(psb_c_ipk_) :: ixb
info = -1;
if (c_associated(x%item)) then
call c_f_pointer(x%item,xp)
else
return
end if
ixb = psb_c_get_index_base()
call xp%set_entry((index+(1-ixb)),val)
info = 0
end function psb_c_zvect_set_entry
function psb_c_zvect_get_entry(x,index) bind(c) result(res)
use psb_base_mod
implicit none
type(psb_c_zvector) :: x
type(psb_z_vect_type), pointer :: xp
integer(psb_c_ipk_), value :: index
complex(c_double_complex) :: res
integer(psb_c_ipk_) :: ixb
if (c_associated(x%item)) then
call c_f_pointer(x%item,xp)
else
return
end if
ixb = psb_c_get_index_base()
res = xp%get_entry((index+(1-ixb)))
end function psb_c_zvect_get_entry
function psb_c_zvect_clone(xh,yh) bind(c) result(info)
implicit none
cbind/base/psb_z_tools_cbind_mod.F90
+11 -5
View File
@@ -281,13 +281,19 @@ contains
end if
ixb = psb_c_get_index_base()
if (ixb == 1) then
call psb_geins(nz,irw(1:nz),val(1:nz),&
& xp,descp,info)
else
select case(ixb)
case (0)
!write(0,*) 'C_GEINS: IDX_BASE',ixb,' :',irw(1:nz)+(1-ixb),val(1:nz)
call psb_geins(nz,(irw(1:nz)+(1-ixb)),val(1:nz),&
& xp,descp,info)
end if
case(1)
!write(0,*) 'C_GEINS: IDX_BASE',ixb,' :',irw(1:nz),val(1:nz)
call psb_geins(nz,irw(1:nz),val(1:nz),&
& xp,descp,info)
case default
write(0,*) 'C_GEINS: Unkonwn inndex base ',ixb
info =-2
end select
res = min(0,info)
cbind/linsolve/psb_base_linsolve_cbind_mod.f90
+15 -1
View File
@@ -15,7 +15,6 @@ contains
implicit none
type(solveroptions) :: options
integer(psb_c_ipk_) :: res
options%itmax = 1000
options%itrace = 0
options%istop = 2
@@ -24,6 +23,21 @@ contains
res = 0
end function psb_c_DefaultSolverOptions
function psb_c_PrintSolverOptions(options)&
& bind(c,name='psb_c_PrintSolverOptions') result(res)
implicit none
type(solveroptions) :: options
integer(psb_c_ipk_) :: res
write(*,*) 'PSBLAS C Interface Solver Options '
write(*,*) ' Maximum number of iterations :', options%itmax
write(*,*) ' Tracing :', options%itrace
write(*,*) ' Stopping Criterion :', options%istop
write(*,*) ' Restart :', options%irst
write(*,*) ' EPS (tolerance) :', options%eps
res = 0
end function psb_c_PrintSolverOptions
end module psb_base_linsolve_cbind_mod
cbind/linsolve/psb_linsolve_cbind.h
+2 -1
View File
@@ -25,7 +25,8 @@ typedef struct psb_c_solveroptions {
} psb_c_SolverOptions;
int psb_c_DefaultSolverOptions(psb_c_SolverOptions *opt);
int psb_c_PrintSolverOptions(psb_c_SolverOptions *opt);
int psb_c_skrylov(const char *method, psb_c_sspmat *ah, psb_c_sprec *ph,
psb_c_svector *bh, psb_c_svector *xh,
psb_c_descriptor *cdh, psb_c_SolverOptions *opt);
config/pac.m4
+2 -20
View File
@@ -2223,6 +2223,7 @@ AS_HELP_STRING([--enable-cuda],
pac_cv_cuda="$enableval";
]
)
AC_MSG_RESULT($pac_cv_cuda)
]
)
@@ -2374,28 +2375,9 @@ AS_HELP_STRING([--enable-openacc],
[Specify whether to enable openacc. ]),
[
pac_cv_openacc="$enableval";
AC_MSG_RESULT([$enableval.])
]
dnl ,
dnl [pac_cv_openacc="no";]
)
if test x"$pac_cv_openacc" == x"yes" ; then
AC_MSG_RESULT([yes.])
# AC_LANG_PUSH([Fortran])
# AC_OPENACC()
# pac_cv_openacc_fcopt="$OPENACC_FCFLAGS";
# AC_LANG_POP()
# AC_LANG_PUSH([C])
# AC_OPENACC()
# pac_cv_openacc_ccopt="$OPENACC_CFLAGS";
# AC_LANG_POP()
# AC_LANG_PUSH([C++])
# AC_OPENACC()
# pac_cv_openacc_cxxopt="$OPENACC_CXXFLAGS";
# AC_LANG_POP()
else
pac_cv_openacc="no";
AC_MSG_RESULT([no.])
fi
]
)
configure
Vendored
+9 -25
View File
@@ -7871,8 +7871,6 @@ fi
if test x"$pac_cv_lpk_size" == x"" ; then
pac_cv_lpk_size=8
fi
PSB_IPKDEF="#define PSB_IPK$pac_cv_ipk_size"
PSB_LPKDEF="#define PSB_LPK$pac_cv_lpk_size"
# Enforce sensible combination
if (( $pac_cv_lpk_size < $pac_cv_ipk_size )); then
{ printf "%s\n" "$as_me:${as_lineno-$LINENO}: Invalid combination of size specs PSB_IPK ${pac_cv_ipk_size} PSB_LPK ${pac_cv_lpk_size}. " >&5
@@ -7881,6 +7879,8 @@ printf "%s\n" "$as_me: Invalid combination of size specs PSB_IPK ${pac_cv_ipk_si
printf "%s\n" "$as_me: Forcing equal values" >&6;}
pac_cv_lpk_size=$pac_cv_ipk_size;
fi
PSB_IPKDEF="#define PSB_IPK$pac_cv_ipk_size"
PSB_LPKDEF="#define PSB_LPK$pac_cv_lpk_size"
FDEFINES="$psblas_cv_define_prepend-DPSB_IPK${pac_cv_ipk_size} $FDEFINES";
FDEFINES="$psblas_cv_define_prepend-DPSB_LPK${pac_cv_lpk_size} $FDEFINES";
@@ -9122,6 +9122,7 @@ else
AR="$ac_cv_prog_AR"
fi
AR="${AR} -cr"
if test -n "$ac_tool_prefix"; then
# Extract the first word of "${ac_tool_prefix}ranlib", so it can be a program name with args.
set dummy ${ac_tool_prefix}ranlib; ac_word=$2
@@ -9226,7 +9227,6 @@ else
RANLIB="$ac_cv_prog_RANLIB"
fi
AR="$AR -cr"
###############################################################################
# BLAS library presence checks
@@ -11287,6 +11287,8 @@ pac_cv_cuda="$enableval";
fi
{ printf "%s\n" "$as_me:${as_lineno-$LINENO}: result: $pac_cv_cuda" >&5
printf "%s\n" "$pac_cv_cuda" >&6; }
if test "x$pac_cv_cuda" == "xyes"; then
@@ -11609,30 +11611,12 @@ if test ${enable_openacc+y}
then :
enableval=$enable_openacc;
pac_cv_openacc="$enableval";
{ printf "%s\n" "$as_me:${as_lineno-$LINENO}: result: $enableval." >&5
printf "%s\n" "$enableval." >&6; }
fi
if test x"$pac_cv_openacc" == x"yes" ; then
{ printf "%s\n" "$as_me:${as_lineno-$LINENO}: result: yes." >&5
printf "%s\n" "yes." >&6; }
# AC_LANG_PUSH([Fortran])
# AC_OPENACC()
# pac_cv_openacc_fcopt="$OPENACC_FCFLAGS";
# AC_LANG_POP()
# AC_LANG_PUSH([C])
# AC_OPENACC()
# pac_cv_openacc_ccopt="$OPENACC_CFLAGS";
# AC_LANG_POP()
# AC_LANG_PUSH([C++])
# AC_OPENACC()
# pac_cv_openacc_cxxopt="$OPENACC_CXXFLAGS";
# AC_LANG_POP()
else
pac_cv_openacc="no";
{ printf "%s\n" "$as_me:${as_lineno-$LINENO}: result: no." >&5
printf "%s\n" "no." >&6; }
fi
if test x"$pac_cv_openacc" == x"yes" ; then
@@ -13523,7 +13507,7 @@ fi
EXTRA_OPT : ${EXTRA_OPT}
CUDA : ${PSB_HAVE_CUDA}
CUDA : ${pac_cv_cuda}
CUDA_CC : ${pac_cv_cudacc}
OPENACC : ${pac_cv_openacc}
@@ -13565,7 +13549,7 @@ printf "%s\n" "$as_me:
EXTRA_OPT : ${EXTRA_OPT}
CUDA : ${PSB_HAVE_CUDA}
CUDA : ${pac_cv_cuda}
CUDA_CC : ${pac_cv_cudacc}
OPENACC : ${pac_cv_openacc}
configure.ac
+5 -5
View File
@@ -563,14 +563,14 @@ fi
if test x"$pac_cv_lpk_size" == x"" ; then
pac_cv_lpk_size=8
fi
PSB_IPKDEF="#define PSB_IPK$pac_cv_ipk_size"
PSB_LPKDEF="#define PSB_LPK$pac_cv_lpk_size"
# Enforce sensible combination
if (( $pac_cv_lpk_size < $pac_cv_ipk_size )); then
AC_MSG_NOTICE([[Invalid combination of size specs PSB_IPK ${pac_cv_ipk_size} PSB_LPK ${pac_cv_lpk_size}. ]]);
AC_MSG_NOTICE([[Forcing equal values]])
pac_cv_lpk_size=$pac_cv_ipk_size;
fi
PSB_IPKDEF="#define PSB_IPK$pac_cv_ipk_size"
PSB_LPKDEF="#define PSB_LPK$pac_cv_lpk_size"
FDEFINES="$psblas_cv_define_prepend-DPSB_IPK${pac_cv_ipk_size} $FDEFINES";
FDEFINES="$psblas_cv_define_prepend-DPSB_LPK${pac_cv_lpk_size} $FDEFINES";
dnl CDEFINES="-DPSB_IPK${pac_cv_ipk_size} -DPSB_LPK${pac_cv_lpk_size} $CDEFINES"
@@ -955,11 +955,11 @@ if test x"$pac_cv_openacc" == x"yes" ; then
dnl CXXOPENACC="$ax_cv_prog_cxx_openacc";
dnl FCOPENACC="$ax_cv_prog_fc_openacc";
dnl else
dnl AC_MSG_NOTICE([OpenACC 1 flags CC $CCOPENACC CXX $CXXOPENACC FC $FCOPENACC])
dnl AC_MSG_NOTICE([OpenACC 1 flags CC $CCOPENACC CXX $CXXOPENACC FC $FCOPENACC])
PAC_ARG_WITH_FLAGS(ccopenacc,CCOPENACC)
PAC_ARG_WITH_FLAGS(cxxopenacc,CXXOPENACC)
PAC_ARG_WITH_FLAGS(fcopenacc,FCOPENACC)
dnl AC_MSG_NOTICE([OpenACC 2 flags CC $CCOPENACC CXX $CXXOPENACC FC $FCOPENACC])
dnl AC_MSG_NOTICE([OpenACC 2 flags CC $CCOPENACC CXX $CXXOPENACC FC $FCOPENACC])
dnl CCOPENACC="$ax_cv_prog_c_openacc";
dnl CXXOPENACC="$ax_cv_prog_cxx_openacc";
dnl FCOPENACC="$ax_cv_prog_fc_openacc";
@@ -1193,7 +1193,7 @@ AC_MSG_NOTICE([
EXTRA_OPT : ${EXTRA_OPT}
CUDA : ${PSB_HAVE_CUDA}
CUDA : ${pac_cv_cuda}
CUDA_CC : ${pac_cv_cudacc}
OPENACC : ${pac_cv_openacc}
cuda/Makefile
+2 -1
View File
@@ -18,7 +18,8 @@ LIBNAME=libpsb_cuda.a
FOBJS=cusparse_mod.o base_cusparse_mod.o \
s_cusparse_mod.o d_cusparse_mod.o c_cusparse_mod.o z_cusparse_mod.o \
s_cusparse_mod.o d_cusparse_mod.o \
c_cusparse_mod.o z_cusparse_mod.o \
psb_vectordev_mod.o core_mod.o \
psb_s_vectordev_mod.o psb_d_vectordev_mod.o psb_i_vectordev_mod.o\
psb_c_vectordev_mod.o psb_z_vectordev_mod.o psb_base_vectordev_mod.o \
cuda/cuda_util.c
+1
View File
@@ -338,6 +338,7 @@ int getGPUMemoryClockRate()
#endif
return(count);
}
int getGPUWarpSize()
{ int count=0;
if (prop!=NULL)
ext/impl/psb_c_hll_csmv.f90
+62 -5
View File
@@ -224,6 +224,24 @@ subroutine psb_c_hll_csmv(alpha,a,x,beta,y,info,trans)
end do
if (info /= psb_success_) goto 9999
case(64)
!$omp parallel do private(i, j,ir,mxrwl, hkpnt)
do i=1,mmhk,hksz
j = ((i-1)/hksz)+1
ir = hksz
mxrwl = (a%hkoffs(j+1) - a%hkoffs(j))/hksz
if (mxrwl>0) then
hkpnt = a%hkoffs(j) + 1
if (info == psb_success_) &
& call psb_c_hll_csmv_notra_64(i,mxrwl,a%irn(i),&
& alpha,a%ja(hkpnt),hksz,a%val(hkpnt),hksz,&
& a%is_triangle(),a%is_unit(),&
& x,beta,y,info)
end if
j = j + 1
end do
if (info /= psb_success_) goto 9999
case default
!$omp parallel do private(i, j,ir,mxrwl, hkpnt)
do i=1,mmhk,hksz
@@ -382,7 +400,7 @@ contains
integer(psb_ipk_), parameter :: m=8
integer(psb_ipk_) :: i,j,k, m4, jc
complex(psb_spk_) :: acc(4), tmp(m)
complex(psb_spk_) :: tmp(m)
info = psb_success_
@@ -420,7 +438,7 @@ contains
integer(psb_ipk_), parameter :: m=24
integer(psb_ipk_) :: i,j,k, m4, jc
complex(psb_spk_) :: acc(4), tmp(m)
complex(psb_spk_) :: tmp(m)
info = psb_success_
@@ -458,7 +476,7 @@ contains
integer(psb_ipk_), parameter :: m=16
integer(psb_ipk_) :: i,j,k, m4, jc
complex(psb_spk_) :: acc(4), tmp(m)
complex(psb_spk_) :: tmp(m)
info = psb_success_
@@ -496,7 +514,7 @@ contains
integer(psb_ipk_), parameter :: m=32
integer(psb_ipk_) :: i,j,k, m4, jc
complex(psb_spk_) :: acc(4), tmp(m)
complex(psb_spk_) :: tmp(m)
info = psb_success_
@@ -522,6 +540,45 @@ contains
end subroutine psb_c_hll_csmv_notra_32
subroutine psb_c_hll_csmv_notra_64(ir,n,irn,alpha,ja,ldj,val,ldv,&
& is_triangle,is_unit, x,beta,y,info)
use psb_base_mod, only : psb_ipk_, psb_spk_, czero, psb_success_
implicit none
integer(psb_ipk_), intent(in) :: ir,n,ldj,ldv,ja(ldj,*),irn(*)
complex(psb_spk_), intent(in) :: alpha, beta, x(*),val(ldv,*)
complex(psb_spk_), intent(inout) :: y(*)
logical, intent(in) :: is_triangle,is_unit
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), parameter :: m=64
integer(psb_ipk_) :: i,j,k, m4, jc
complex(psb_spk_) :: tmp(m)
info = psb_success_
tmp(:) = czero
if (alpha /= czero) then
do j=1, maxval(irn(1:64))
tmp(1:64) = tmp(1:64) + val(1:64,j)*x(ja(1:64,j))
end do
end if
if (beta == czero) then
y(ir:ir+64-1) = alpha*tmp(1:64)
else
y(ir:ir+64-1) = alpha*tmp(1:64) + beta*y(ir:ir+64-1)
end if
if (is_unit) then
do i=1, min(64,n)
y(ir+i-1) = y(ir+i-1) + alpha*x(ir+i-1)
end do
end if
end subroutine psb_c_hll_csmv_notra_64
subroutine psb_c_hll_csmv_notra_4(ir,n,irn,alpha,ja,ldj,val,ldv,&
& is_triangle,is_unit, x,beta,y,info)
use psb_base_mod, only : psb_ipk_, psb_spk_, czero, psb_success_
@@ -534,7 +591,7 @@ contains
integer(psb_ipk_), parameter :: m=4
integer(psb_ipk_) :: i,j,k, m4, jc
complex(psb_spk_) :: acc(4), tmp(m)
complex(psb_spk_) :: tmp(m)
info = psb_success_
ext/impl/psb_d_hll_csmv.f90
+62 -5
View File
@@ -224,6 +224,24 @@ subroutine psb_d_hll_csmv(alpha,a,x,beta,y,info,trans)
end do
if (info /= psb_success_) goto 9999
case(64)
!$omp parallel do private(i, j,ir,mxrwl, hkpnt)
do i=1,mmhk,hksz
j = ((i-1)/hksz)+1
ir = hksz
mxrwl = (a%hkoffs(j+1) - a%hkoffs(j))/hksz
if (mxrwl>0) then
hkpnt = a%hkoffs(j) + 1
if (info == psb_success_) &
& call psb_d_hll_csmv_notra_64(i,mxrwl,a%irn(i),&
& alpha,a%ja(hkpnt),hksz,a%val(hkpnt),hksz,&
& a%is_triangle(),a%is_unit(),&
& x,beta,y,info)
end if
j = j + 1
end do
if (info /= psb_success_) goto 9999
case default
!$omp parallel do private(i, j,ir,mxrwl, hkpnt)
do i=1,mmhk,hksz
@@ -382,7 +400,7 @@ contains
integer(psb_ipk_), parameter :: m=8
integer(psb_ipk_) :: i,j,k, m4, jc
real(psb_dpk_) :: acc(4), tmp(m)
real(psb_dpk_) :: tmp(m)
info = psb_success_
@@ -420,7 +438,7 @@ contains
integer(psb_ipk_), parameter :: m=24
integer(psb_ipk_) :: i,j,k, m4, jc
real(psb_dpk_) :: acc(4), tmp(m)
real(psb_dpk_) :: tmp(m)
info = psb_success_
@@ -458,7 +476,7 @@ contains
integer(psb_ipk_), parameter :: m=16
integer(psb_ipk_) :: i,j,k, m4, jc
real(psb_dpk_) :: acc(4), tmp(m)
real(psb_dpk_) :: tmp(m)
info = psb_success_
@@ -496,7 +514,7 @@ contains
integer(psb_ipk_), parameter :: m=32
integer(psb_ipk_) :: i,j,k, m4, jc
real(psb_dpk_) :: acc(4), tmp(m)
real(psb_dpk_) :: tmp(m)
info = psb_success_
@@ -522,6 +540,45 @@ contains
end subroutine psb_d_hll_csmv_notra_32
subroutine psb_d_hll_csmv_notra_64(ir,n,irn,alpha,ja,ldj,val,ldv,&
& is_triangle,is_unit, x,beta,y,info)
use psb_base_mod, only : psb_ipk_, psb_dpk_, dzero, psb_success_
implicit none
integer(psb_ipk_), intent(in) :: ir,n,ldj,ldv,ja(ldj,*),irn(*)
real(psb_dpk_), intent(in) :: alpha, beta, x(*),val(ldv,*)
real(psb_dpk_), intent(inout) :: y(*)
logical, intent(in) :: is_triangle,is_unit
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), parameter :: m=64
integer(psb_ipk_) :: i,j,k, m4, jc
real(psb_dpk_) :: tmp(m)
info = psb_success_
tmp(:) = dzero
if (alpha /= dzero) then
do j=1, maxval(irn(1:64))
tmp(1:64) = tmp(1:64) + val(1:64,j)*x(ja(1:64,j))
end do
end if
if (beta == dzero) then
y(ir:ir+64-1) = alpha*tmp(1:64)
else
y(ir:ir+64-1) = alpha*tmp(1:64) + beta*y(ir:ir+64-1)
end if
if (is_unit) then
do i=1, min(64,n)
y(ir+i-1) = y(ir+i-1) + alpha*x(ir+i-1)
end do
end if
end subroutine psb_d_hll_csmv_notra_64
subroutine psb_d_hll_csmv_notra_4(ir,n,irn,alpha,ja,ldj,val,ldv,&
& is_triangle,is_unit, x,beta,y,info)
use psb_base_mod, only : psb_ipk_, psb_dpk_, dzero, psb_success_
@@ -534,7 +591,7 @@ contains
integer(psb_ipk_), parameter :: m=4
integer(psb_ipk_) :: i,j,k, m4, jc
real(psb_dpk_) :: acc(4), tmp(m)
real(psb_dpk_) :: tmp(m)
info = psb_success_
ext/impl/psb_s_hll_csmv.f90
+62 -5
View File
@@ -224,6 +224,24 @@ subroutine psb_s_hll_csmv(alpha,a,x,beta,y,info,trans)
end do
if (info /= psb_success_) goto 9999
case(64)
!$omp parallel do private(i, j,ir,mxrwl, hkpnt)
do i=1,mmhk,hksz
j = ((i-1)/hksz)+1
ir = hksz
mxrwl = (a%hkoffs(j+1) - a%hkoffs(j))/hksz
if (mxrwl>0) then
hkpnt = a%hkoffs(j) + 1
if (info == psb_success_) &
& call psb_s_hll_csmv_notra_64(i,mxrwl,a%irn(i),&
& alpha,a%ja(hkpnt),hksz,a%val(hkpnt),hksz,&
& a%is_triangle(),a%is_unit(),&
& x,beta,y,info)
end if
j = j + 1
end do
if (info /= psb_success_) goto 9999
case default
!$omp parallel do private(i, j,ir,mxrwl, hkpnt)
do i=1,mmhk,hksz
@@ -382,7 +400,7 @@ contains
integer(psb_ipk_), parameter :: m=8
integer(psb_ipk_) :: i,j,k, m4, jc
real(psb_spk_) :: acc(4), tmp(m)
real(psb_spk_) :: tmp(m)
info = psb_success_
@@ -420,7 +438,7 @@ contains
integer(psb_ipk_), parameter :: m=24
integer(psb_ipk_) :: i,j,k, m4, jc
real(psb_spk_) :: acc(4), tmp(m)
real(psb_spk_) :: tmp(m)
info = psb_success_
@@ -458,7 +476,7 @@ contains
integer(psb_ipk_), parameter :: m=16
integer(psb_ipk_) :: i,j,k, m4, jc
real(psb_spk_) :: acc(4), tmp(m)
real(psb_spk_) :: tmp(m)
info = psb_success_
@@ -496,7 +514,7 @@ contains
integer(psb_ipk_), parameter :: m=32
integer(psb_ipk_) :: i,j,k, m4, jc
real(psb_spk_) :: acc(4), tmp(m)
real(psb_spk_) :: tmp(m)
info = psb_success_
@@ -522,6 +540,45 @@ contains
end subroutine psb_s_hll_csmv_notra_32
subroutine psb_s_hll_csmv_notra_64(ir,n,irn,alpha,ja,ldj,val,ldv,&
& is_triangle,is_unit, x,beta,y,info)
use psb_base_mod, only : psb_ipk_, psb_spk_, szero, psb_success_
implicit none
integer(psb_ipk_), intent(in) :: ir,n,ldj,ldv,ja(ldj,*),irn(*)
real(psb_spk_), intent(in) :: alpha, beta, x(*),val(ldv,*)
real(psb_spk_), intent(inout) :: y(*)
logical, intent(in) :: is_triangle,is_unit
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), parameter :: m=64
integer(psb_ipk_) :: i,j,k, m4, jc
real(psb_spk_) :: tmp(m)
info = psb_success_
tmp(:) = szero
if (alpha /= szero) then
do j=1, maxval(irn(1:64))
tmp(1:64) = tmp(1:64) + val(1:64,j)*x(ja(1:64,j))
end do
end if
if (beta == szero) then
y(ir:ir+64-1) = alpha*tmp(1:64)
else
y(ir:ir+64-1) = alpha*tmp(1:64) + beta*y(ir:ir+64-1)
end if
if (is_unit) then
do i=1, min(64,n)
y(ir+i-1) = y(ir+i-1) + alpha*x(ir+i-1)
end do
end if
end subroutine psb_s_hll_csmv_notra_64
subroutine psb_s_hll_csmv_notra_4(ir,n,irn,alpha,ja,ldj,val,ldv,&
& is_triangle,is_unit, x,beta,y,info)
use psb_base_mod, only : psb_ipk_, psb_spk_, szero, psb_success_
@@ -534,7 +591,7 @@ contains
integer(psb_ipk_), parameter :: m=4
integer(psb_ipk_) :: i,j,k, m4, jc
real(psb_spk_) :: acc(4), tmp(m)
real(psb_spk_) :: tmp(m)
info = psb_success_
ext/impl/psb_z_hll_csmv.f90
+62 -5
View File
@@ -224,6 +224,24 @@ subroutine psb_z_hll_csmv(alpha,a,x,beta,y,info,trans)
end do
if (info /= psb_success_) goto 9999
case(64)
!$omp parallel do private(i, j,ir,mxrwl, hkpnt)
do i=1,mmhk,hksz
j = ((i-1)/hksz)+1
ir = hksz
mxrwl = (a%hkoffs(j+1) - a%hkoffs(j))/hksz
if (mxrwl>0) then
hkpnt = a%hkoffs(j) + 1
if (info == psb_success_) &
& call psb_z_hll_csmv_notra_64(i,mxrwl,a%irn(i),&
& alpha,a%ja(hkpnt),hksz,a%val(hkpnt),hksz,&
& a%is_triangle(),a%is_unit(),&
& x,beta,y,info)
end if
j = j + 1
end do
if (info /= psb_success_) goto 9999
case default
!$omp parallel do private(i, j,ir,mxrwl, hkpnt)
do i=1,mmhk,hksz
@@ -382,7 +400,7 @@ contains
integer(psb_ipk_), parameter :: m=8
integer(psb_ipk_) :: i,j,k, m4, jc
complex(psb_dpk_) :: acc(4), tmp(m)
complex(psb_dpk_) :: tmp(m)
info = psb_success_
@@ -420,7 +438,7 @@ contains
integer(psb_ipk_), parameter :: m=24
integer(psb_ipk_) :: i,j,k, m4, jc
complex(psb_dpk_) :: acc(4), tmp(m)
complex(psb_dpk_) :: tmp(m)
info = psb_success_
@@ -458,7 +476,7 @@ contains
integer(psb_ipk_), parameter :: m=16
integer(psb_ipk_) :: i,j,k, m4, jc
complex(psb_dpk_) :: acc(4), tmp(m)
complex(psb_dpk_) :: tmp(m)
info = psb_success_
@@ -496,7 +514,7 @@ contains
integer(psb_ipk_), parameter :: m=32
integer(psb_ipk_) :: i,j,k, m4, jc
complex(psb_dpk_) :: acc(4), tmp(m)
complex(psb_dpk_) :: tmp(m)
info = psb_success_
@@ -522,6 +540,45 @@ contains
end subroutine psb_z_hll_csmv_notra_32
subroutine psb_z_hll_csmv_notra_64(ir,n,irn,alpha,ja,ldj,val,ldv,&
& is_triangle,is_unit, x,beta,y,info)
use psb_base_mod, only : psb_ipk_, psb_dpk_, zzero, psb_success_
implicit none
integer(psb_ipk_), intent(in) :: ir,n,ldj,ldv,ja(ldj,*),irn(*)
complex(psb_dpk_), intent(in) :: alpha, beta, x(*),val(ldv,*)
complex(psb_dpk_), intent(inout) :: y(*)
logical, intent(in) :: is_triangle,is_unit
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), parameter :: m=64
integer(psb_ipk_) :: i,j,k, m4, jc
complex(psb_dpk_) :: tmp(m)
info = psb_success_
tmp(:) = zzero
if (alpha /= zzero) then
do j=1, maxval(irn(1:64))
tmp(1:64) = tmp(1:64) + val(1:64,j)*x(ja(1:64,j))
end do
end if
if (beta == zzero) then
y(ir:ir+64-1) = alpha*tmp(1:64)
else
y(ir:ir+64-1) = alpha*tmp(1:64) + beta*y(ir:ir+64-1)
end if
if (is_unit) then
do i=1, min(64,n)
y(ir+i-1) = y(ir+i-1) + alpha*x(ir+i-1)
end do
end if
end subroutine psb_z_hll_csmv_notra_64
subroutine psb_z_hll_csmv_notra_4(ir,n,irn,alpha,ja,ldj,val,ldv,&
& is_triangle,is_unit, x,beta,y,info)
use psb_base_mod, only : psb_ipk_, psb_dpk_, zzero, psb_success_
@@ -534,7 +591,7 @@ contains
integer(psb_ipk_), parameter :: m=4
integer(psb_ipk_) :: i,j,k, m4, jc
complex(psb_dpk_) :: acc(4), tmp(m)
complex(psb_dpk_) :: tmp(m)
info = psb_success_
linsolve/impl/CMakeLists.txt
+4
View File
@@ -42,6 +42,10 @@ set(PSB_krylov_source_files
psb_z_krylov_conv_mod.f90
psb_zkrylov.f90
psb_zrgmres.f90
psb_dminres.f90
psb_sminres.f90
psb_cminres.f90
psb_zminres.f90
)
foreach(file IN LISTS PSB_krylov_source_files)
list(APPEND krylov_source_files ${CMAKE_CURRENT_LIST_DIR}/${file})
linsolve/impl/Makefile
+2 -1
View File
@@ -13,7 +13,8 @@ OBJS=psb_dkrylov.o psb_skrylov.o psb_ckrylov.o psb_zkrylov.o \
psb_ccgstab.o psb_ccg.o psb_cfcg.o psb_cgcr.o psb_ccgs.o \
psb_cbicg.o psb_ccgstabl.o psb_crgmres.o\
psb_zcgstab.o psb_zcg.o psb_zfcg.o psb_zgcr.o psb_zcgs.o \
psb_zbicg.o psb_zcgstabl.o psb_zrgmres.o
psb_zbicg.o psb_zcgstabl.o psb_zrgmres.o \
psb_dminres.o psb_sminres.o psb_cminres.o psb_zminres.o
LIBNAME=$(METHDLIBNAME)
COBJS=
linsolve/impl/psb_cbicg.f90
+18 -8
View File
@@ -160,19 +160,29 @@ subroutine psb_cbicg_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
endif
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if(info /= psb_success_) then
linsolve/impl/psb_ccg.F90
+22 -1
View File
@@ -159,8 +159,29 @@ subroutine psb_ccg_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_)&
linsolve/impl/psb_ccgs.f90
+22 -1
View File
@@ -154,8 +154,29 @@ Subroutine psb_ccgs_vect(a,prec,b,x,eps,desc_a,info,&
If (Present(istop)) Then
istop_ = istop
Else
istop_ = 2
istop_ = psb_get_istop_default()
Endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_) call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info)
linsolve/impl/psb_ccgstab.f90
+23 -5
View File
@@ -156,13 +156,31 @@ Subroutine psb_ccgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
If (Present(istop)) Then
istop_ = istop
Else
istop_ = 2
Endif
else
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! ISTOP_ = 1: Normwise backward error, infinity norm
! ISTOP_ = 2: ||r||/||b|| norm 2
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
! = if (.not.same_type_as(x,b)) then
! = write(0,*) 'Warning: different dynamic types for X and B '
! = end if
linsolve/impl/psb_ccgstabl.f90
+22 -1
View File
@@ -172,8 +172,29 @@ Subroutine psb_ccgstabl_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
if (present(itmax)) then
itmax_ = itmax
linsolve/impl/psb_cfcg.F90
+22 -2
View File
@@ -164,9 +164,29 @@ subroutine psb_cfcg_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_)&
linsolve/impl/psb_cgcr.f90
+19 -11
View File
@@ -167,22 +167,30 @@ subroutine psb_cgcr_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
!
! ISTOP_ = 1: Normwise backward error, infinity norm
! ISTOP_ = 2: ||r||/||b||, 2-norm
!
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
endif
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_)&
& call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info)
linsolve/impl/psb_ckrylov.f90
+4 -1
View File
@@ -150,7 +150,7 @@ Subroutine psb_ckrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
procedure(psb_ckryl_vect) :: psb_cbicg_vect, psb_ccgstab_vect,&
& psb_ccgs_vect
procedure(psb_ckryl_rest_vect) :: psb_crgmres_vect, psb_ccgstabl_vect, psb_cgcr_vect
procedure(psb_ckryl_cond_vect) :: psb_ccg_vect, psb_cfcg_vect
procedure(psb_ckryl_cond_vect) :: psb_ccg_vect, psb_cfcg_vect, psb_cminres_vect
logical :: do_alloc_wrk
type(psb_ctxt_type) :: ctxt
@@ -199,6 +199,9 @@ Subroutine psb_ckrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
case('RGMRES','GMRES')
call psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,irst=irst,istop=istop)
case('MINRES','PMINRES')
call psb_cminres_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop)
case('BICGSTABL')
call psb_ccgstabl_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,irst=irst,istop=istop)
linsolve/impl/psb_cminres.f90
+511
View File
@@ -0,0 +1,511 @@
!
! Parallel Sparse BLAS version 3.5
! (C) Copyright 2006-2018
! Salvatore Filippone
! Alfredo Buttari
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
!
! File: psb_cminres.f90
! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
! C C
! C References: C
! C [1] Duff, I., Marrone, M., Radicati, G., and Vittoli, C. C
! C Level 3 basic linear algebra subprograms for sparse C
! C matrices: a user level interface C
! C ACM Trans. Math. Softw., 23(3), 379-401, 1997. C
! C C
! C C
! C [2] S. Filippone, M. Colajanni C
! C PSBLAS: A library for parallel linear algebra C
! C computation on sparse matrices C
! C ACM Trans. on Math. Softw., 26(4), 527-550, Dec. 2000. C
! C C
! C [3] M. Arioli, I. Duff, M. Ruiz C
! C Stopping criteria for iterative solvers C
! C SIAM J. Matrix Anal. Appl., Vol. 13, pp. 138-144, 1992 C
! C C
! C C
! C [4] R. Barrett et al C
! C Templates for the solution of linear systems C
! C SIAM, 1993
! C C
! C C
! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
! File: psb_cminres.f90
!
! Subroutine: psb_cminres
! This subroutine implements the MINRES method.
!
!
! Arguments:
!
! a - type(psb_cspmat_type) Input: sparse matrix containing A.
! prec - class(psb_cprec_type) Input: preconditioner
! b(:) - complex Input: vector containing the
! right hand side B
! x(:) - complex Input/Output: vector containing the
! initial guess and final solution X.
! eps - real Input: Stopping tolerance; the iteration is
! stopped when the error estimate |err| <= eps
! desc_a - type(psb_desc_type). Input: The communication descriptor.
! info - integer. Output: Return code
!
! itmax - integer(optional) Input: maximum number of iterations to be
! performed.
! iter - integer(optional) Output: how many iterations have been
! performed.
! performed.
! err - real (optional) Output: error estimate on exit. If the
! denominator of the estimate is exactly
! 0, it is changed into 1.
! itrace - integer(optional) Input: print an informational message
! with the error estimate every itrace
! iterations
! istop - integer(optional) Input: stopping criterion, or how
! to estimate the error.
! 1: err = |r|/(|a||x|+|b|); here the iteration is
! stopped when |r| <= eps * (|a||x|+|b|)
! 2: err = |r|/|b|; here the iteration is
! stopped when |r| <= eps * |b|
! where r is the (preconditioned, recursive
! estimate of) residual.
!
!
subroutine psb_cminres_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop,cond)
use psb_base_mod
use psb_prec_mod
use psb_c_linsolve_conv_mod
use psb_linsolve_mod
implicit none
type(psb_cspmat_type), intent(in) :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(psb_cprec_type), intent(inout) :: prec
type(psb_c_vect_type), Intent(inout) :: b
type(psb_c_vect_type), Intent(inout) :: x
Real(psb_spk_), Intent(in) :: eps
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_spk_), Optional, Intent(out) :: err, cond
! = Local data
complex(psb_spk_), allocatable, target :: aux(:)
type(psb_c_vect_type), allocatable, target :: wwrk(:)
type(psb_c_vect_type), pointer :: y, r1, r2, v, w, w1, w2, wtmp, res
complex(psb_spk_) :: alfa, beta, beta1, dbar, delta, denom
complex(psb_spk_) :: epsln, gamma, gbar, oldb, oldeps, phi, phibar
complex(psb_spk_) :: rhs1, rhs2, s, zeta, betatmp
complex(psb_spk_) :: rotx(1), roty(1)
real(psb_spk_) :: cs
complex(psb_spk_) :: sn
integer(psb_ipk_) :: itmax_, istop_, naux, itx, itrace_, n_col, n_row, err_act
integer(psb_lpk_) :: mglob
integer(psb_ipk_) :: debug_level, debug_unit
type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: np, me
real(psb_dpk_) :: derr, errnum, errden, deps
logical :: do_cond
logical :: converged
real(psb_spk_) :: epsmach
real(psb_spk_) :: gmax, gmin, tnorm2, ynorm2
real(psb_spk_) :: rni, xni, bni, ani, bn2, r0n2
character(len=20) :: name
character(len=*), parameter :: methdname='MINRES'
info = psb_success_
name = 'psb_cminres'
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ctxt = desc_a%get_context()
call psb_info(ctxt, me, np)
if (.not.allocated(b%v)) then
info = psb_err_invalid_vect_state_
call psb_errpush(info,name)
goto 9999
endif
if (.not.allocated(x%v)) then
info = psb_err_invalid_vect_state_
call psb_errpush(info,name)
goto 9999
endif
mglob = desc_a%get_global_rows()
n_row = desc_a%get_local_rows()
n_col = desc_a%get_local_cols()
if (present(istop)) then
istop_ = istop
else
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_)&
& call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_chkvect on X/B')
goto 9999
end if
naux=4*n_col
allocate(aux(naux), stat=info)
if (info == psb_success_) call psb_geall(wwrk,desc_a,info,n=9_psb_ipk_)
if (info == psb_success_) call psb_geasb(wwrk,desc_a,info,mold=x%v,scratch=.true.)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
if (present(itmax)) then
itmax_ = itmax
else
itmax_ = 1000
endif
if (present(itrace)) then
itrace_ = itrace
else
itrace_ = 0
end if
y => wwrk(1)
r1 => wwrk(2)
r2 => wwrk(3)
v => wwrk(4)
w => wwrk(5)
w1 => wwrk(6)
w2 => wwrk(7)
wtmp => wwrk(8)
res => wwrk(9)
do_cond = present(cond)
epsmach = epsilon(sone)
! res = b - A*x
call psb_geaxpby(cone,b,czero,res,desc_a,info)
if (info == psb_success_) call psb_spmm(-cone,a,x,cone,res,desc_a,info,work=aux)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
select case(istop_)
case(psb_istop_ani_)
ani = psb_spnrmi(a,desc_a,info)
bni = psb_geamax(b,desc_a,info)
case(psb_istop_bn2_)
bn2 = psb_genrm2(b,desc_a,info)
case(psb_istop_rn2_abs_)
! nothing needed
case(psb_istop_rrn2_)
r0n2 = psb_genrm2(res,desc_a,info)
end select
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
errnum = czero
errden = cone
deps = eps
converged = .false.
if ((itrace_ > 0).and.(me == 0)) call log_header(methdname)
! y = beta1 * P' * v1, with v1 the first Lanczos vector.
call psb_geaxpby(cone,res,czero,y,desc_a,info)
if (info == psb_success_) call psb_geaxpby(cone,res,czero,r1,desc_a,info)
if (info == psb_success_) call prec%apply(res,y,desc_a,info,work=aux)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
beta1 = sqrt(psb_gedot(res,y,desc_a,info))
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
! If beta1 is numerically zero, x is already acceptable.
if (abs(beta1) <= epsmach) then
itx = 0
converged = .true.
if (do_cond) cond = szero
select case(istop_)
case(psb_istop_ani_)
rni = psb_geamax(res,desc_a,info)
xni = psb_geamax(x,desc_a,info)
errnum = rni
errden = ani*xni + bni
case(psb_istop_bn2_)
errnum = abs(beta1)
errden = bn2
case(psb_istop_rn2_abs_)
errnum = abs(beta1)
errden = sone
case(psb_istop_rrn2_)
errnum = abs(beta1)
errden = r0n2
end select
if (errden == szero) errden = sone
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
end if
beta = beta1
oldb = czero
dbar = czero
epsln = czero
phibar = beta1
rhs1 = beta1
rhs2 = czero
cs = -sone
sn = czero
itx = 0
if (do_cond) then
tnorm2 = 0.0_psb_spk_
ynorm2 = 0.0_psb_spk_
gmax = 0.0_psb_spk_
gmin = huge(sone)
end if
call psb_geaxpby(cone,r1,czero,r2,desc_a,info)
if (info == psb_success_) call psb_geaxpby(czero,x,czero,w,desc_a,info)
if (info == psb_success_) call psb_geaxpby(czero,x,czero,w1,desc_a,info)
if (info == psb_success_) call psb_geaxpby(czero,x,czero,w2,desc_a,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
if (.not.converged) then
do while (itx < itmax_)
if (abs(beta) <= epsmach) exit
itx = itx + 1
s = cone/beta
call psb_geaxpby(s,y,czero,v,desc_a,info)
if (info == psb_success_) call psb_spmm(cone,a,v,czero,y,desc_a,info,work=aux)
if (itx >= 2 .and. info == psb_success_) then
call psb_geaxpby((-beta/oldb),r1,cone,y,desc_a,info)
end if
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
alfa = psb_gedot(v,y,desc_a,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
call psb_geaxpby((-alfa/beta),r2,cone,y,desc_a,info)
if (info == psb_success_) call psb_geaxpby(cone,r2,czero,r1,desc_a,info)
if (info == psb_success_) call psb_geaxpby(cone,y,czero,r2,desc_a,info)
if (info == psb_success_) call prec%apply(r2,y,desc_a,info,work=aux)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
oldb = beta
beta = sqrt(psb_gedot(r2,y,desc_a,info))
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
if (do_cond) then
tnorm2 = tnorm2 + abs(alfa)**2 + abs(oldb)**2 + abs(beta)**2
if (itx == 1) then
gmax = abs(alfa)
gmin = gmax
end if
end if
oldeps = epsln
delta = dbar
gbar = alfa
rotx(1) = delta
roty(1) = gbar
call crot(1,rotx,1,roty,1,cs,sn)
delta = rotx(1)
gbar = -roty(1)
epsln = czero
dbar = beta
rotx(1) = epsln
roty(1) = dbar
call crot(1,rotx,1,roty,1,cs,sn)
epsln = rotx(1)
dbar = -roty(1)
gamma = gbar
betatmp = beta
call crotg(gamma,betatmp,cs,sn)
if (abs(gamma) <= epsmach) exit
phi = phibar
phibar = czero
rotx(1) = phi
roty(1) = phibar
call crot(1,rotx,1,roty,1,cs,-sn)
phi = rotx(1)
phibar = roty(1)
denom = cone/gamma
call psb_geaxpby(cone,w2,czero,w1,desc_a,info)
if (info == psb_success_) call psb_geaxpby(cone,w,czero,w2,desc_a,info)
if (info == psb_success_) call psb_geaxpby(cone,v,czero,wtmp,desc_a,info)
if (info == psb_success_) call psb_geaxpby(-oldeps,w1,cone,wtmp,desc_a,info)
if (info == psb_success_) call psb_geaxpby(-delta,w2,cone,wtmp,desc_a,info)
if (info == psb_success_) call psb_geaxpby(denom,wtmp,czero,w,desc_a,info)
if (info == psb_success_) call psb_geaxpby(phi,w,cone,x,desc_a,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
if (do_cond) then
gmax = max(gmax,abs(gamma))
gmin = min(gmin,abs(gamma))
zeta = rhs1/gamma
ynorm2 = ynorm2 + abs(zeta)**2
rhs1 = rhs2 - delta*zeta
rhs2 = -epsln*zeta
end if
select case(istop_)
case(psb_istop_ani_)
! Compute true residual only for the ANI stopping criterion.
call psb_geaxpby(cone,b,czero,res,desc_a,info)
if (info == psb_success_) call psb_spmm(-cone,a,x,cone,res,desc_a,info,work=aux)
if (info == psb_success_) rni = psb_geamax(res,desc_a,info)
if (info == psb_success_) xni = psb_geamax(x,desc_a,info)
errnum = rni
errden = ani*xni + bni
case(psb_istop_bn2_)
errnum = abs(phibar)
errden = bn2
case(psb_istop_rn2_abs_)
errnum = abs(phibar)
errden = sone
case(psb_istop_rrn2_)
errnum = abs(phibar)
errden = r0n2
end select
if (errden == szero) errden = sone
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,deps)
if (errnum <= eps*errden) then
converged = .true.
exit
end if
end do
end if
if ((.not.converged) .and. (itx >= itmax_)) then
if (debug_level >= psb_debug_ext_)&
& write(debug_unit,*) me,' ',trim(name),': MINRES reached iteration limit'
end if
if (do_cond) then
if (gmin > 0.0_psb_spk_) then
cond = gmax/gmin
else
cond = huge(sone)
end if
end if
call log_end(methdname,me,itx,itrace_,errnum,errden,deps,err=derr,iter=iter)
if (present(err)) err = derr
if (info == psb_success_) call psb_gefree(wwrk,desc_a,info)
if (info == psb_success_) deallocate(aux,stat=info)
if (info /= psb_success_) then
call psb_errpush(info,name)
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine psb_cminres_vect
linsolve/impl/psb_crgmres.f90
+133 -89
View File
@@ -97,7 +97,8 @@
! iterations
! istop - integer(optional) Input: stopping criterion, or how
! to estimate the error.
! 1: err = |r|/(|a||x|+|b|); here the iteration is
! 1: err = |r|/(|a||x|+|b|); here
! the iteration is
! stopped when |r| <= eps * (|a||x|+|b|)
! 2: err = |r|/|b|; here the iteration is
! stopped when |r| <= eps * |b|
@@ -130,7 +131,7 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
type(psb_c_vect_type) :: w, w1, xt
real(psb_spk_) :: tmp
complex(psb_spk_) :: scal, gm, rti, rti1
integer(psb_ipk_) ::litmax, naux, it, k, itrace_,&
integer(psb_ipk_) ::itmax_, naux, it, k, itrace_,&
& n_row, n_col, nl
integer(psb_lpk_) :: mglob
Logical, Parameter :: exchange=.True., noexchange=.False., use_srot=.true.
@@ -172,24 +173,35 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
!
! ISTOP_ = 1: Normwise backward error, infinity norm
! ISTOP_ = 2: ||r||/||b||, 2-norm
!
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
endif
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_,psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
if (present(itmax)) then
litmax = itmax
itmax_ = itmax
else
litmax = 1000
itmax_ = 1000
endif
if (present(itrace)) then
@@ -253,12 +265,15 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
& w%get_nrows(),w1%get_nrows()
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
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)
else if (istop_ == 3) then
case(psb_istop_bn2_)
bn2 = psb_genrm2(b,desc_a,info)
case(psb_istop_rn2_abs_)
! do nothing
case(psb_istop_rrn2_)
call psb_geaxpby(cone,b,czero,v(1),desc_a,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
@@ -273,7 +288,8 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
goto 9999
end if
r0n2 = psb_genrm2(v(1),desc_a,info)
endif
end select
errnum = czero
errden = cone
deps = eps
@@ -324,27 +340,32 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
!
! check convergence
!
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
rni = psb_geamax(v(1),desc_a,info)
xni = psb_geamax(x,desc_a,info)
errnum = rni
errden = (ani*xni+bni)
else if (istop_ == 2) then
case(psb_istop_bn2_)
rni = psb_genrm2(v(1),desc_a,info)
errnum = rni
errden = bn2
else if (istop_ == 3) then
case(psb_istop_rn2_abs_)
rni = psb_genrm2(v(1),desc_a,info)
errnum = rni
errden = sone
case(psb_istop_rrn2_)
rni = psb_genrm2(v(1),desc_a,info)
errnum = rni
errden = r0n2
endif
end select
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
if ((errnum <= eps*errden).or.(itx >= litmax)) exit restart
if ((errnum <= eps*errden).or.(itx >= itmax_)) exit restart
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,deps)
@@ -360,42 +381,18 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
call prec%apply(v(i),w1,desc_a,info)
call psb_spmm(cone,a,w1,czero,w,desc_a,info,work=aux)
!
call mgs(i,h,v,w,rs,c,s,desc_a,info)
do k = 1, i
h(k,i) = psb_gedot(v(k),w,desc_a,info)
call psb_geaxpby(-h(k,i),v(k),cone,w,desc_a,info)
end do
h(i+1,i) = psb_genrm2(w,desc_a,info)
scal=cone/h(i+1,i)
call psb_geaxpby(scal,w,czero,v(i+1),desc_a,info)
do k=2,i
call crot(1,h(k-1,i),1,h(k,i),1,real(c(k-1),kind=psb_spk_),s(k-1))
enddo
rti = h(i,i)
rti1 = h(i+1,i)
call crotg(rti,rti1,tmp,s(i))
c(i) = cmplx(tmp,szero)
call crot(1,h(i,i),1,h(i+1,i),1,real(c(i),kind=psb_spk_),s(i))
h(i+1,i) = czero
call crot(1,rs(i),1,rs(i+1),1,real(c(i),kind=psb_spk_),s(i))
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
!
! build x and then compute the residual and its infinity norm
!
rst = rs
call w1%set(czero)
call ctrsm('l','u','n','n',i,1,cone,h,size(h,1),rst,size(rst,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rst(1:i)
do k=1, i
call psb_geaxpby(rst(k),v(k),cone,xt,desc_a,info)
end do
call prec%apply(xt,desc_a,info)
call psb_geaxpby(cone,x,cone,xt,desc_a,info)
call psb_geaxpby(cone,x,czero,xt,desc_a,info)
call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info)
call psb_geaxpby(cone,b,czero,w1,desc_a,info)
call psb_spmm(-cone,a,xt,cone,w1,desc_a,info,work=aux)
rni = psb_geamax(w1,desc_a,info)
@@ -403,8 +400,7 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
errnum = rni
errden = (ani*xni+bni)
!
else if (istop_ == 2) then
case(psb_istop_bn2_)
!
! compute the residual 2-norm as byproduct of the solution
! procedure of the least-squares problem
@@ -412,7 +408,14 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
rni = abs(rs(i+1))
errnum = rni
errden = bn2
else if (istop_ == 3) then
case(psb_istop_rn2_abs_)
rni = abs(rs(i+1))
errnum = rni
errden = sone
case(psb_istop_rrn2_)
!
! compute the residual 2-norm as byproduct of the solution
! procedure of the least-squares problem
@@ -420,28 +423,18 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
rni = abs(rs(i+1))
errnum = rni
errden = r0n2
endif
end select
if (errnum <= eps*errden) then
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
call psb_geaxpby(cone,xt,czero,x,desc_a,info)
! = x = xt
else if (istop_ == 2) then
!
! build x
!
call ctrsm('l','u','n','n',i,1,cone,h,size(h,1),rs,size(rs,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:i)
call w1%set(czero)
do k=1, i
call psb_geaxpby(rs(k),v(k),cone,w1,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(cone,w,cone,x,desc_a,info)
end if
case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_)
call rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info) !
end select
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,ione,errnum,errden,deps)
@@ -454,25 +447,15 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
end do inner
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
call psb_geaxpby(cone,xt,czero,x,desc_a,info)! x = xt
else if (istop_ == 2) then
!
! build x
!
call ctrsm('l','u','n','n',nl,1,cone,h,size(h,1),rs,size(rs,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:nl)
call w1%set(czero)
do k=1, nl
call psb_geaxpby(rs(k),v(k),cone,w1,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(cone,w,cone,x,desc_a,info)
end if
case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_)
call rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info) !
end select
if (itx >= litmax) then
if (itx >= itmax_) then
if (itrace_ > 0) then
if (mod(itx,itrace_)/=0) &
& call log_conv(methdname,me,itx,ione,errnum,errden,deps)
@@ -502,6 +485,67 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
9999 call psb_error_handler(err_act)
return
contains
!
! Advance one step the Gram-Schmidt process, and apply
! Givens rotations to keep H triangular
!
subroutine mgs(n,h,v,w,rs,c,s,desc_a,info)
complex(psb_spk_) :: c(:), s(:), h(:,:), rs(:)
type(psb_c_vect_type) :: v(:), w
type(psb_desc_type) :: desc_a
complex(psb_spk_) :: scal, gm, rti, rti1
real(psb_spk_) :: tmp
integer(psb_ipk_) :: info
integer(psb_ipk_) :: k,n
!
do k = 1, n
h(k,n) = psb_gedot(v(k),w,desc_a,info)
call psb_geaxpby(-h(k,n),v(k),cone,w,desc_a,info)
end do
h(n+1,n) = psb_genrm2(w,desc_a,info)
scal=cone/h(n+1,n)
call psb_geaxpby(scal,w,czero,v(n+1),desc_a,info)
do k=2,n
call crot(1,h(k-1:k-1,n),1,h(k:k,n),1,real(c(k-1),kind=psb_spk_),s(k-1))
enddo
rti = h(n,n)
rti1 = h(n+1,n)
call crotg(rti,rti1,tmp,s(n))
c(n) = cmplx(tmp,szero)
call crot(1,h(n:n,n),1,h(n+1:n+1,n),1,real(c(n),kind=psb_spk_),s(n))
call crot(1,rs(n:n),1,rs(n+1:n+1),1,real(c(n),kind=psb_spk_),s(n))
end subroutine mgs
!
! Rebuild solution X from the space V using the factor
! stored in R
!
subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info)
complex(psb_spk_) :: c(:), s(:), rs(:), h(:,:)
type(psb_c_vect_type) :: v(:), w, w1, x
type(psb_desc_type) :: desc_a
class(psb_cprec_type) :: prec
integer(psb_ipk_) :: info
integer(psb_ipk_) :: k,n
!
!
! build x
!
call ctrsm('l','u','n','n',n,1,cone,h,size(h,1),rs,size(rs,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:n)
call w1%zero()
do k=1, n
call psb_geaxpby(rs(k),v(k),cone,w1,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(cone,w,cone,x,desc_a,info)
end subroutine rebuildx
end subroutine psb_crgmres_vect
linsolve/impl/psb_dbicg.f90
+18 -8
View File
@@ -160,19 +160,29 @@ subroutine psb_dbicg_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
endif
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if(info /= psb_success_) then
linsolve/impl/psb_dcg.F90
+22 -1
View File
@@ -159,8 +159,29 @@ subroutine psb_dcg_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_)&
linsolve/impl/psb_dcgs.f90
+22 -1
View File
@@ -154,8 +154,29 @@ Subroutine psb_dcgs_vect(a,prec,b,x,eps,desc_a,info,&
If (Present(istop)) Then
istop_ = istop
Else
istop_ = 2
istop_ = psb_get_istop_default()
Endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_) call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info)
linsolve/impl/psb_dcgstab.f90
+23 -5
View File
@@ -156,13 +156,31 @@ Subroutine psb_dcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
If (Present(istop)) Then
istop_ = istop
Else
istop_ = 2
Endif
else
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! ISTOP_ = 1: Normwise backward error, infinity norm
! ISTOP_ = 2: ||r||/||b|| norm 2
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
! = if (.not.same_type_as(x,b)) then
! = write(0,*) 'Warning: different dynamic types for X and B '
! = end if
linsolve/impl/psb_dcgstabl.f90
+22 -1
View File
@@ -172,8 +172,29 @@ Subroutine psb_dcgstabl_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
if (present(itmax)) then
itmax_ = itmax
linsolve/impl/psb_dfcg.F90
+22 -2
View File
@@ -164,9 +164,29 @@ subroutine psb_dfcg_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_)&
linsolve/impl/psb_dgcr.f90
+19 -11
View File
@@ -167,22 +167,30 @@ subroutine psb_dgcr_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
!
! ISTOP_ = 1: Normwise backward error, infinity norm
! ISTOP_ = 2: ||r||/||b||, 2-norm
!
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
endif
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_)&
& call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info)
linsolve/impl/psb_dkrylov.f90
+4 -1
View File
@@ -150,7 +150,7 @@ Subroutine psb_dkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
procedure(psb_dkryl_vect) :: psb_dbicg_vect, psb_dcgstab_vect,&
& psb_dcgs_vect
procedure(psb_dkryl_rest_vect) :: psb_drgmres_vect, psb_dcgstabl_vect, psb_dgcr_vect
procedure(psb_dkryl_cond_vect) :: psb_dcg_vect, psb_dfcg_vect
procedure(psb_dkryl_cond_vect) :: psb_dcg_vect, psb_dfcg_vect, psb_dminres_vect
logical :: do_alloc_wrk
type(psb_ctxt_type) :: ctxt
@@ -199,6 +199,9 @@ Subroutine psb_dkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
case('RGMRES','GMRES')
call psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,irst=irst,istop=istop)
case('MINRES','PMINRES')
call psb_dminres_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop)
case('BICGSTABL')
call psb_dcgstabl_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,irst=irst,istop=istop)
linsolve/impl/psb_dminres.f90
+526
View File
@@ -0,0 +1,526 @@
!
! Parallel Sparse BLAS version 3.5
! (C) Copyright 2006-2018
! Salvatore Filippone
! Alfredo Buttari
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
!
! File: psb_dminres.f90
! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
! C C
! C References: C
! C [1] Duff, I., Marrone, M., Radicati, G., and Vittoli, C. C
! C Level 3 basic linear algebra subprograms for sparse C
! C matrices: a user level interface C
! C ACM Trans. Math. Softw., 23(3), 379-401, 1997. C
! C C
! C C
! C [2] S. Filippone, M. Colajanni C
! C PSBLAS: A library for parallel linear algebra C
! C computation on sparse matrices C
! C ACM Trans. on Math. Softw., 26(4), 527-550, Dec. 2000. C
! C C
! C [3] M. Arioli, I. Duff, M. Ruiz C
! C Stopping criteria for iterative solvers C
! C SIAM J. Matrix Anal. Appl., Vol. 13, pp. 138-144, 1992 C
! C C
! C C
! C [4] R. Barrett et al C
! C Templates for the solution of linear systems C
! C SIAM, 1993
! C C
! C C
! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
! File: psb_dminres.f90
!
! Subroutine: psb_dminres
! This subroutine implements the MINRES method.
!
!
! Arguments:
!
! a - type(psb_dspmat_type) Input: sparse matrix containing A.
! prec - class(psb_dprec_type) Input: preconditioner
! b(:) - real Input: vector containing the
! right hand side B
! x(:) - real Input/Output: vector containing the
! initial guess and final solution X.
! eps - real Input: Stopping tolerance; the iteration is
! stopped when the error estimate |err| <= eps
! desc_a - type(psb_desc_type). Input: The communication descriptor.
! info - integer. Output: Return code
!
! itmax - integer(optional) Input: maximum number of iterations to be
! performed.
! iter - integer(optional) Output: how many iterations have been
! performed.
! performed.
! err - real (optional) Output: error estimate on exit. If the
! denominator of the estimate is exactly
! 0, it is changed into 1.
! itrace - integer(optional) Input: print an informational message
! with the error estimate every itrace
! iterations
! istop - integer(optional) Input: stopping criterion, or how
! to estimate the error.
! 1: err = |r|/(|a||x|+|b|); here the iteration is
! stopped when |r| <= eps * (|a||x|+|b|)
! 2: err = |r|/|b|; here the iteration is
! stopped when |r| <= eps * |b|
! where r is the (preconditioned, recursive
! estimate of) residual.
!
!
subroutine psb_dminres_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop,cond)
use psb_base_mod
use psb_prec_mod
use psb_d_linsolve_conv_mod
use psb_linsolve_mod
implicit none
type(psb_dspmat_type), intent(in) :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(psb_dprec_type), intent(inout) :: prec
type(psb_d_vect_type), Intent(inout) :: b
type(psb_d_vect_type), Intent(inout) :: x
Real(psb_dpk_), Intent(in) :: eps
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_dpk_), Optional, Intent(out) :: err, cond
! = Local data
real(psb_dpk_), allocatable, target :: aux(:)
type(psb_d_vect_type), allocatable, target :: wwrk(:)
type(psb_d_vect_type), pointer :: y, r1, r2, v, w, w1, w2, wtmp, res
real(psb_dpk_) :: alfa, beta, beta1, dbar, delta, denom
real(psb_dpk_) :: epsln, gamma, gbar, oldb, oldeps, phi, phibar
real(psb_dpk_) :: rhs1, rhs2, s, zeta, betatmp
real(psb_dpk_) :: rotx(1), roty(1)
real(psb_dpk_) :: cs
real(psb_dpk_) :: sn
integer(psb_ipk_) :: itmax_, istop_, naux, itx, itrace_, n_col, n_row, err_act
integer(psb_lpk_) :: mglob
integer(psb_ipk_) :: debug_level, debug_unit
type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: np, me
real(psb_dpk_) :: derr, errnum, errden, deps
logical :: do_cond
logical :: converged
real(psb_dpk_) :: epsmach
real(psb_dpk_) :: gmax, gmin, tnorm2, ynorm2
real(psb_dpk_) :: rni, xni, bni, ani, bn2, r0n2
real(psb_dpk_) :: dotv
character(len=20) :: name
character(len=*), parameter :: methdname='MINRES'
info = psb_success_
name = 'psb_dminres'
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ctxt = desc_a%get_context()
call psb_info(ctxt, me, np)
if (.not.allocated(b%v)) then
info = psb_err_invalid_vect_state_
call psb_errpush(info,name)
goto 9999
endif
if (.not.allocated(x%v)) then
info = psb_err_invalid_vect_state_
call psb_errpush(info,name)
goto 9999
endif
mglob = desc_a%get_global_rows()
n_row = desc_a%get_local_rows()
n_col = desc_a%get_local_cols()
if (present(istop)) then
istop_ = istop
else
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_)&
& call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_chkvect on X/B')
goto 9999
end if
naux=4*n_col
allocate(aux(naux), stat=info)
if (info == psb_success_) call psb_geall(wwrk,desc_a,info,n=9_psb_ipk_)
if (info == psb_success_) call psb_geasb(wwrk,desc_a,info,mold=x%v,scratch=.true.)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
if (present(itmax)) then
itmax_ = itmax
else
itmax_ = 1000
endif
if (present(itrace)) then
itrace_ = itrace
else
itrace_ = 0
end if
y => wwrk(1)
r1 => wwrk(2)
r2 => wwrk(3)
v => wwrk(4)
w => wwrk(5)
w1 => wwrk(6)
w2 => wwrk(7)
wtmp => wwrk(8)
res => wwrk(9)
do_cond = present(cond)
epsmach = epsilon(done)
! res = b - A*x
call psb_geaxpby(done,b,dzero,res,desc_a,info)
if (info == psb_success_) call psb_spmm(-done,a,x,done,res,desc_a,info,work=aux)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
select case(istop_)
case(psb_istop_ani_)
ani = psb_spnrmi(a,desc_a,info)
bni = psb_geamax(b,desc_a,info)
case(psb_istop_bn2_)
bn2 = psb_genrm2(b,desc_a,info)
case(psb_istop_rn2_abs_)
! nothing needed
case(psb_istop_rrn2_)
r0n2 = psb_genrm2(res,desc_a,info)
end select
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
errnum = dzero
errden = done
deps = eps
converged = .false.
if ((itrace_ > 0).and.(me == 0)) call log_header(methdname)
! y = beta1 * P' * v1, with v1 the first Lanczos vector.
call psb_geaxpby(done,res,dzero,y,desc_a,info)
if (info == psb_success_) call psb_geaxpby(done,res,dzero,r1,desc_a,info)
if (info == psb_success_) call prec%apply(res,y,desc_a,info,work=aux)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
dotv = psb_gedot(res,y,desc_a,info)
if (dotv < -epsmach*max(done,abs(dotv))) then
info=psb_err_pivot_too_small_
call psb_errpush(info,name,a_err='MINRES breakdown: negative (r,z)')
goto 9999
end if
if (dotv < dzero) dotv = dzero
beta1 = sqrt(dotv)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
! If beta1 is numerically zero, x is already acceptable.
if (abs(beta1) <= epsmach) then
itx = 0
converged = .true.
if (do_cond) cond = dzero
select case(istop_)
case(psb_istop_ani_)
rni = psb_geamax(res,desc_a,info)
xni = psb_geamax(x,desc_a,info)
errnum = rni
errden = ani*xni + bni
case(psb_istop_bn2_)
errnum = abs(beta1)
errden = bn2
case(psb_istop_rn2_abs_)
errnum = abs(beta1)
errden = done
case(psb_istop_rrn2_)
errnum = abs(beta1)
errden = r0n2
end select
if (errden == dzero) errden = done
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
end if
beta = beta1
oldb = dzero
dbar = dzero
epsln = dzero
phibar = beta1
rhs1 = beta1
rhs2 = dzero
cs = -done
sn = dzero
itx = 0
if (do_cond) then
tnorm2 = 0.0_psb_dpk_
ynorm2 = 0.0_psb_dpk_
gmax = 0.0_psb_dpk_
gmin = huge(done)
end if
call psb_geaxpby(done,r1,dzero,r2,desc_a,info)
if (info == psb_success_) call psb_geaxpby(dzero,x,dzero,w,desc_a,info)
if (info == psb_success_) call psb_geaxpby(dzero,x,dzero,w1,desc_a,info)
if (info == psb_success_) call psb_geaxpby(dzero,x,dzero,w2,desc_a,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
if (.not.converged) then
do while (itx < itmax_)
if (abs(beta) <= epsmach) exit
itx = itx + 1
s = done/beta
call psb_geaxpby(s,y,dzero,v,desc_a,info)
if (info == psb_success_) call psb_spmm(done,a,v,dzero,y,desc_a,info,work=aux)
if (itx >= 2 .and. info == psb_success_) then
call psb_geaxpby((-beta/oldb),r1,done,y,desc_a,info)
end if
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
alfa = psb_gedot(v,y,desc_a,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
call psb_geaxpby((-alfa/beta),r2,done,y,desc_a,info)
if (info == psb_success_) call psb_geaxpby(done,r2,dzero,r1,desc_a,info)
if (info == psb_success_) call psb_geaxpby(done,y,dzero,r2,desc_a,info)
if (info == psb_success_) call prec%apply(r2,y,desc_a,info,work=aux)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
oldb = beta
dotv = psb_gedot(r2,y,desc_a,info)
if (dotv < -epsmach*max(done,abs(dotv))) then
info=psb_err_pivot_too_small_
call psb_errpush(info,name,a_err='MINRES breakdown: negative (r2,z2)')
goto 9999
end if
if (dotv < dzero) dotv = dzero
beta = sqrt(dotv)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
if (do_cond) then
tnorm2 = tnorm2 + abs(alfa)**2 + abs(oldb)**2 + abs(beta)**2
if (itx == 1) then
gmax = abs(alfa)
gmin = gmax
end if
end if
oldeps = epsln
delta = dbar
gbar = alfa
rotx(1) = delta
roty(1) = gbar
call drot(1,rotx,1,roty,1,cs,sn)
delta = rotx(1)
gbar = -roty(1)
epsln = dzero
dbar = beta
rotx(1) = epsln
roty(1) = dbar
call drot(1,rotx,1,roty,1,cs,sn)
epsln = rotx(1)
dbar = -roty(1)
gamma = gbar
betatmp = beta
call drotg(gamma,betatmp,cs,sn)
if (abs(gamma) <= epsmach) exit
phi = phibar
phibar = dzero
rotx(1) = phi
roty(1) = phibar
call drot(1,rotx,1,roty,1,cs,-sn)
phi = rotx(1)
phibar = roty(1)
denom = done/gamma
call psb_geaxpby(done,w2,dzero,w1,desc_a,info)
if (info == psb_success_) call psb_geaxpby(done,w,dzero,w2,desc_a,info)
if (info == psb_success_) call psb_geaxpby(done,v,dzero,wtmp,desc_a,info)
if (info == psb_success_) call psb_geaxpby(-oldeps,w1,done,wtmp,desc_a,info)
if (info == psb_success_) call psb_geaxpby(-delta,w2,done,wtmp,desc_a,info)
if (info == psb_success_) call psb_geaxpby(denom,wtmp,dzero,w,desc_a,info)
if (info == psb_success_) call psb_geaxpby(phi,w,done,x,desc_a,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
if (do_cond) then
gmax = max(gmax,abs(gamma))
gmin = min(gmin,abs(gamma))
zeta = rhs1/gamma
ynorm2 = ynorm2 + abs(zeta)**2
rhs1 = rhs2 - delta*zeta
rhs2 = -epsln*zeta
end if
select case(istop_)
case(psb_istop_ani_)
! Compute true residual only for the ANI stopping criterion.
call psb_geaxpby(done,b,dzero,res,desc_a,info)
if (info == psb_success_) call psb_spmm(-done,a,x,done,res,desc_a,info,work=aux)
if (info == psb_success_) rni = psb_geamax(res,desc_a,info)
if (info == psb_success_) xni = psb_geamax(x,desc_a,info)
errnum = rni
errden = ani*xni + bni
case(psb_istop_bn2_)
errnum = abs(phibar)
errden = bn2
case(psb_istop_rn2_abs_)
errnum = abs(phibar)
errden = done
case(psb_istop_rrn2_)
errnum = abs(phibar)
errden = r0n2
end select
if (errden == dzero) errden = done
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,deps)
if (errnum <= eps*errden) then
converged = .true.
exit
end if
end do
end if
if ((.not.converged) .and. (itx >= itmax_)) then
if (debug_level >= psb_debug_ext_)&
& write(debug_unit,*) me,' ',trim(name),': MINRES reached iteration limit'
end if
if (do_cond) then
if (gmin > 0.0_psb_dpk_) then
cond = gmax/gmin
else
cond = huge(done)
end if
end if
call log_end(methdname,me,itx,itrace_,errnum,errden,deps,err=derr,iter=iter)
if (present(err)) err = derr
if (info == psb_success_) call psb_gefree(wwrk,desc_a,info)
if (info == psb_success_) deallocate(aux,stat=info)
if (info /= psb_success_) then
call psb_errpush(info,name)
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine psb_dminres_vect
linsolve/impl/psb_drgmres.f90
+133 -89
View File
@@ -97,7 +97,8 @@
! iterations
! istop - integer(optional) Input: stopping criterion, or how
! to estimate the error.
! 1: err = |r|/(|a||x|+|b|); here the iteration is
! 1: err = |r|/(|a||x|+|b|); here
! the iteration is
! stopped when |r| <= eps * (|a||x|+|b|)
! 2: err = |r|/|b|; here the iteration is
! stopped when |r| <= eps * |b|
@@ -130,7 +131,7 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
type(psb_d_vect_type) :: w, w1, xt
real(psb_dpk_) :: tmp
real(psb_dpk_) :: scal, gm, rti, rti1
integer(psb_ipk_) ::litmax, naux, it, k, itrace_,&
integer(psb_ipk_) ::itmax_, naux, it, k, itrace_,&
& n_row, n_col, nl
integer(psb_lpk_) :: mglob
Logical, Parameter :: exchange=.True., noexchange=.False., use_srot=.true.
@@ -172,24 +173,35 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
!
! ISTOP_ = 1: Normwise backward error, infinity norm
! ISTOP_ = 2: ||r||/||b||, 2-norm
!
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
endif
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_,psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
if (present(itmax)) then
litmax = itmax
itmax_ = itmax
else
litmax = 1000
itmax_ = 1000
endif
if (present(itrace)) then
@@ -253,12 +265,15 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
& w%get_nrows(),w1%get_nrows()
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
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)
else if (istop_ == 3) then
case(psb_istop_bn2_)
bn2 = psb_genrm2(b,desc_a,info)
case(psb_istop_rn2_abs_)
! do nothing
case(psb_istop_rrn2_)
call psb_geaxpby(done,b,dzero,v(1),desc_a,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
@@ -273,7 +288,8 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
goto 9999
end if
r0n2 = psb_genrm2(v(1),desc_a,info)
endif
end select
errnum = dzero
errden = done
deps = eps
@@ -324,27 +340,32 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
!
! check convergence
!
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
rni = psb_geamax(v(1),desc_a,info)
xni = psb_geamax(x,desc_a,info)
errnum = rni
errden = (ani*xni+bni)
else if (istop_ == 2) then
case(psb_istop_bn2_)
rni = psb_genrm2(v(1),desc_a,info)
errnum = rni
errden = bn2
else if (istop_ == 3) then
case(psb_istop_rn2_abs_)
rni = psb_genrm2(v(1),desc_a,info)
errnum = rni
errden = done
case(psb_istop_rrn2_)
rni = psb_genrm2(v(1),desc_a,info)
errnum = rni
errden = r0n2
endif
end select
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
if ((errnum <= eps*errden).or.(itx >= litmax)) exit restart
if ((errnum <= eps*errden).or.(itx >= itmax_)) exit restart
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,deps)
@@ -360,42 +381,18 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
call prec%apply(v(i),w1,desc_a,info)
call psb_spmm(done,a,w1,dzero,w,desc_a,info,work=aux)
!
call mgs(i,h,v,w,rs,c,s,desc_a,info)
do k = 1, i
h(k,i) = psb_gedot(v(k),w,desc_a,info)
call psb_geaxpby(-h(k,i),v(k),done,w,desc_a,info)
end do
h(i+1,i) = psb_genrm2(w,desc_a,info)
scal=done/h(i+1,i)
call psb_geaxpby(scal,w,dzero,v(i+1),desc_a,info)
do k=2,i
call drot(1,h(k-1,i),1,h(k,i),1,real(c(k-1),kind=psb_dpk_),s(k-1))
enddo
rti = h(i,i)
rti1 = h(i+1,i)
call drotg(rti,rti1,tmp,s(i))
c(i) = cmplx(tmp,dzero)
call drot(1,h(i,i),1,h(i+1,i),1,real(c(i),kind=psb_dpk_),s(i))
h(i+1,i) = dzero
call drot(1,rs(i),1,rs(i+1),1,real(c(i),kind=psb_dpk_),s(i))
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
!
! build x and then compute the residual and its infinity norm
!
rst = rs
call w1%set(dzero)
call dtrsm('l','u','n','n',i,1,done,h,size(h,1),rst,size(rst,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rst(1:i)
do k=1, i
call psb_geaxpby(rst(k),v(k),done,xt,desc_a,info)
end do
call prec%apply(xt,desc_a,info)
call psb_geaxpby(done,x,done,xt,desc_a,info)
call psb_geaxpby(done,x,dzero,xt,desc_a,info)
call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info)
call psb_geaxpby(done,b,dzero,w1,desc_a,info)
call psb_spmm(-done,a,xt,done,w1,desc_a,info,work=aux)
rni = psb_geamax(w1,desc_a,info)
@@ -403,8 +400,7 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
errnum = rni
errden = (ani*xni+bni)
!
else if (istop_ == 2) then
case(psb_istop_bn2_)
!
! compute the residual 2-norm as byproduct of the solution
! procedure of the least-squares problem
@@ -412,7 +408,14 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
rni = abs(rs(i+1))
errnum = rni
errden = bn2
else if (istop_ == 3) then
case(psb_istop_rn2_abs_)
rni = abs(rs(i+1))
errnum = rni
errden = done
case(psb_istop_rrn2_)
!
! compute the residual 2-norm as byproduct of the solution
! procedure of the least-squares problem
@@ -420,28 +423,18 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
rni = abs(rs(i+1))
errnum = rni
errden = r0n2
endif
end select
if (errnum <= eps*errden) then
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
call psb_geaxpby(done,xt,dzero,x,desc_a,info)
! = x = xt
else if (istop_ == 2) then
!
! build x
!
call dtrsm('l','u','n','n',i,1,done,h,size(h,1),rs,size(rs,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:i)
call w1%set(dzero)
do k=1, i
call psb_geaxpby(rs(k),v(k),done,w1,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(done,w,done,x,desc_a,info)
end if
case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_)
call rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info) !
end select
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,ione,errnum,errden,deps)
@@ -454,25 +447,15 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
end do inner
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
call psb_geaxpby(done,xt,dzero,x,desc_a,info)! x = xt
else if (istop_ == 2) then
!
! build x
!
call dtrsm('l','u','n','n',nl,1,done,h,size(h,1),rs,size(rs,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:nl)
call w1%set(dzero)
do k=1, nl
call psb_geaxpby(rs(k),v(k),done,w1,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(done,w,done,x,desc_a,info)
end if
case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_)
call rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info) !
end select
if (itx >= litmax) then
if (itx >= itmax_) then
if (itrace_ > 0) then
if (mod(itx,itrace_)/=0) &
& call log_conv(methdname,me,itx,ione,errnum,errden,deps)
@@ -502,6 +485,67 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
9999 call psb_error_handler(err_act)
return
contains
!
! Advance one step the Gram-Schmidt process, and apply
! Givens rotations to keep H triangular
!
subroutine mgs(n,h,v,w,rs,c,s,desc_a,info)
real(psb_dpk_) :: c(:), s(:), h(:,:), rs(:)
type(psb_d_vect_type) :: v(:), w
type(psb_desc_type) :: desc_a
real(psb_dpk_) :: scal, gm, rti, rti1
real(psb_dpk_) :: tmp
integer(psb_ipk_) :: info
integer(psb_ipk_) :: k,n
!
do k = 1, n
h(k,n) = psb_gedot(v(k),w,desc_a,info)
call psb_geaxpby(-h(k,n),v(k),done,w,desc_a,info)
end do
h(n+1,n) = psb_genrm2(w,desc_a,info)
scal=done/h(n+1,n)
call psb_geaxpby(scal,w,dzero,v(n+1),desc_a,info)
do k=2,n
call drot(1,h(k-1:k-1,n),1,h(k:k,n),1,real(c(k-1),kind=psb_dpk_),s(k-1))
enddo
rti = h(n,n)
rti1 = h(n+1,n)
call drotg(rti,rti1,tmp,s(n))
c(n) = cmplx(tmp,dzero)
call drot(1,h(n:n,n),1,h(n+1:n+1,n),1,real(c(n),kind=psb_dpk_),s(n))
call drot(1,rs(n:n),1,rs(n+1:n+1),1,real(c(n),kind=psb_dpk_),s(n))
end subroutine mgs
!
! Rebuild solution X from the space V using the factor
! stored in R
!
subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info)
real(psb_dpk_) :: c(:), s(:), rs(:), h(:,:)
type(psb_d_vect_type) :: v(:), w, w1, x
type(psb_desc_type) :: desc_a
class(psb_dprec_type) :: prec
integer(psb_ipk_) :: info
integer(psb_ipk_) :: k,n
!
!
! build x
!
call dtrsm('l','u','n','n',n,1,done,h,size(h,1),rs,size(rs,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:n)
call w1%zero()
do k=1, n
call psb_geaxpby(rs(k),v(k),done,w1,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(done,w,done,x,desc_a,info)
end subroutine rebuildx
end subroutine psb_drgmres_vect
linsolve/impl/psb_sbicg.f90
+18 -8
View File
@@ -160,19 +160,29 @@ subroutine psb_sbicg_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
endif
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if(info /= psb_success_) then
linsolve/impl/psb_scg.F90
+22 -1
View File
@@ -159,8 +159,29 @@ subroutine psb_scg_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_)&
linsolve/impl/psb_scgs.f90
+22 -1
View File
@@ -154,8 +154,29 @@ Subroutine psb_scgs_vect(a,prec,b,x,eps,desc_a,info,&
If (Present(istop)) Then
istop_ = istop
Else
istop_ = 2
istop_ = psb_get_istop_default()
Endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_) call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info)
linsolve/impl/psb_scgstab.f90
+23 -5
View File
@@ -156,13 +156,31 @@ Subroutine psb_scgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
If (Present(istop)) Then
istop_ = istop
Else
istop_ = 2
Endif
else
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! ISTOP_ = 1: Normwise backward error, infinity norm
! ISTOP_ = 2: ||r||/||b|| norm 2
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
! = if (.not.same_type_as(x,b)) then
! = write(0,*) 'Warning: different dynamic types for X and B '
! = end if
linsolve/impl/psb_scgstabl.f90
+22 -1
View File
@@ -172,8 +172,29 @@ Subroutine psb_scgstabl_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
if (present(itmax)) then
itmax_ = itmax
linsolve/impl/psb_sfcg.F90
+22 -2
View File
@@ -164,9 +164,29 @@ subroutine psb_sfcg_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_)&
linsolve/impl/psb_sgcr.f90
+19 -11
View File
@@ -167,22 +167,30 @@ subroutine psb_sgcr_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
!
! ISTOP_ = 1: Normwise backward error, infinity norm
! ISTOP_ = 2: ||r||/||b||, 2-norm
!
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
endif
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_)&
& call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info)
linsolve/impl/psb_skrylov.f90
+4 -1
View File
@@ -150,7 +150,7 @@ Subroutine psb_skrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
procedure(psb_skryl_vect) :: psb_sbicg_vect, psb_scgstab_vect,&
& psb_scgs_vect
procedure(psb_skryl_rest_vect) :: psb_srgmres_vect, psb_scgstabl_vect, psb_sgcr_vect
procedure(psb_skryl_cond_vect) :: psb_scg_vect, psb_sfcg_vect
procedure(psb_skryl_cond_vect) :: psb_scg_vect, psb_sfcg_vect, psb_sminres_vect
logical :: do_alloc_wrk
type(psb_ctxt_type) :: ctxt
@@ -199,6 +199,9 @@ Subroutine psb_skrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
case('RGMRES','GMRES')
call psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,irst=irst,istop=istop)
case('MINRES','PMINRES')
call psb_sminres_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop)
case('BICGSTABL')
call psb_scgstabl_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,irst=irst,istop=istop)
linsolve/impl/psb_sminres.f90
+526
View File
@@ -0,0 +1,526 @@
!
! Parallel Sparse BLAS version 3.5
! (C) Copyright 2006-2018
! Salvatore Filippone
! Alfredo Buttari
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
!
! File: psb_sminres.f90
! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
! C C
! C References: C
! C [1] Duff, I., Marrone, M., Radicati, G., and Vittoli, C. C
! C Level 3 basic linear algebra subprograms for sparse C
! C matrices: a user level interface C
! C ACM Trans. Math. Softw., 23(3), 379-401, 1997. C
! C C
! C C
! C [2] S. Filippone, M. Colajanni C
! C PSBLAS: A library for parallel linear algebra C
! C computation on sparse matrices C
! C ACM Trans. on Math. Softw., 26(4), 527-550, Dec. 2000. C
! C C
! C [3] M. Arioli, I. Duff, M. Ruiz C
! C Stopping criteria for iterative solvers C
! C SIAM J. Matrix Anal. Appl., Vol. 13, pp. 138-144, 1992 C
! C C
! C C
! C [4] R. Barrett et al C
! C Templates for the solution of linear systems C
! C SIAM, 1993
! C C
! C C
! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
! File: psb_sminres.f90
!
! Subroutine: psb_sminres
! This subroutine implements the MINRES method.
!
!
! Arguments:
!
! a - type(psb_sspmat_type) Input: sparse matrix containing A.
! prec - class(psb_sprec_type) Input: preconditioner
! b(:) - real Input: vector containing the
! right hand side B
! x(:) - real Input/Output: vector containing the
! initial guess and final solution X.
! eps - real Input: Stopping tolerance; the iteration is
! stopped when the error estimate |err| <= eps
! desc_a - type(psb_desc_type). Input: The communication descriptor.
! info - integer. Output: Return code
!
! itmax - integer(optional) Input: maximum number of iterations to be
! performed.
! iter - integer(optional) Output: how many iterations have been
! performed.
! performed.
! err - real (optional) Output: error estimate on exit. If the
! denominator of the estimate is exactly
! 0, it is changed into 1.
! itrace - integer(optional) Input: print an informational message
! with the error estimate every itrace
! iterations
! istop - integer(optional) Input: stopping criterion, or how
! to estimate the error.
! 1: err = |r|/(|a||x|+|b|); here the iteration is
! stopped when |r| <= eps * (|a||x|+|b|)
! 2: err = |r|/|b|; here the iteration is
! stopped when |r| <= eps * |b|
! where r is the (preconditioned, recursive
! estimate of) residual.
!
!
subroutine psb_sminres_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop,cond)
use psb_base_mod
use psb_prec_mod
use psb_s_linsolve_conv_mod
use psb_linsolve_mod
implicit none
type(psb_sspmat_type), intent(in) :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(psb_sprec_type), intent(inout) :: prec
type(psb_s_vect_type), Intent(inout) :: b
type(psb_s_vect_type), Intent(inout) :: x
Real(psb_spk_), Intent(in) :: eps
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_spk_), Optional, Intent(out) :: err, cond
! = Local data
real(psb_spk_), allocatable, target :: aux(:)
type(psb_s_vect_type), allocatable, target :: wwrk(:)
type(psb_s_vect_type), pointer :: y, r1, r2, v, w, w1, w2, wtmp, res
real(psb_spk_) :: alfa, beta, beta1, dbar, delta, denom
real(psb_spk_) :: epsln, gamma, gbar, oldb, oldeps, phi, phibar
real(psb_spk_) :: rhs1, rhs2, s, zeta, betatmp
real(psb_spk_) :: rotx(1), roty(1)
real(psb_spk_) :: cs
real(psb_spk_) :: sn
integer(psb_ipk_) :: itmax_, istop_, naux, itx, itrace_, n_col, n_row, err_act
integer(psb_lpk_) :: mglob
integer(psb_ipk_) :: debug_level, debug_unit
type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: np, me
real(psb_dpk_) :: derr, errnum, errden, deps
logical :: do_cond
logical :: converged
real(psb_spk_) :: epsmach
real(psb_spk_) :: gmax, gmin, tnorm2, ynorm2
real(psb_spk_) :: rni, xni, bni, ani, bn2, r0n2
real(psb_spk_) :: dotv
character(len=20) :: name
character(len=*), parameter :: methdname='MINRES'
info = psb_success_
name = 'psb_sminres'
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ctxt = desc_a%get_context()
call psb_info(ctxt, me, np)
if (.not.allocated(b%v)) then
info = psb_err_invalid_vect_state_
call psb_errpush(info,name)
goto 9999
endif
if (.not.allocated(x%v)) then
info = psb_err_invalid_vect_state_
call psb_errpush(info,name)
goto 9999
endif
mglob = desc_a%get_global_rows()
n_row = desc_a%get_local_rows()
n_col = desc_a%get_local_cols()
if (present(istop)) then
istop_ = istop
else
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_)&
& call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_chkvect on X/B')
goto 9999
end if
naux=4*n_col
allocate(aux(naux), stat=info)
if (info == psb_success_) call psb_geall(wwrk,desc_a,info,n=9_psb_ipk_)
if (info == psb_success_) call psb_geasb(wwrk,desc_a,info,mold=x%v,scratch=.true.)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
if (present(itmax)) then
itmax_ = itmax
else
itmax_ = 1000
endif
if (present(itrace)) then
itrace_ = itrace
else
itrace_ = 0
end if
y => wwrk(1)
r1 => wwrk(2)
r2 => wwrk(3)
v => wwrk(4)
w => wwrk(5)
w1 => wwrk(6)
w2 => wwrk(7)
wtmp => wwrk(8)
res => wwrk(9)
do_cond = present(cond)
epsmach = epsilon(sone)
! res = b - A*x
call psb_geaxpby(sone,b,szero,res,desc_a,info)
if (info == psb_success_) call psb_spmm(-sone,a,x,sone,res,desc_a,info,work=aux)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
select case(istop_)
case(psb_istop_ani_)
ani = psb_spnrmi(a,desc_a,info)
bni = psb_geamax(b,desc_a,info)
case(psb_istop_bn2_)
bn2 = psb_genrm2(b,desc_a,info)
case(psb_istop_rn2_abs_)
! nothing needed
case(psb_istop_rrn2_)
r0n2 = psb_genrm2(res,desc_a,info)
end select
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
errnum = szero
errden = sone
deps = eps
converged = .false.
if ((itrace_ > 0).and.(me == 0)) call log_header(methdname)
! y = beta1 * P' * v1, with v1 the first Lanczos vector.
call psb_geaxpby(sone,res,szero,y,desc_a,info)
if (info == psb_success_) call psb_geaxpby(sone,res,szero,r1,desc_a,info)
if (info == psb_success_) call prec%apply(res,y,desc_a,info,work=aux)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
dotv = psb_gedot(res,y,desc_a,info)
if (dotv < -epsmach*max(sone,abs(dotv))) then
info=psb_err_pivot_too_small_
call psb_errpush(info,name,a_err='MINRES breakdown: negative (r,z)')
goto 9999
end if
if (dotv < szero) dotv = szero
beta1 = sqrt(dotv)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
! If beta1 is numerically zero, x is already acceptable.
if (abs(beta1) <= epsmach) then
itx = 0
converged = .true.
if (do_cond) cond = szero
select case(istop_)
case(psb_istop_ani_)
rni = psb_geamax(res,desc_a,info)
xni = psb_geamax(x,desc_a,info)
errnum = rni
errden = ani*xni + bni
case(psb_istop_bn2_)
errnum = abs(beta1)
errden = bn2
case(psb_istop_rn2_abs_)
errnum = abs(beta1)
errden = sone
case(psb_istop_rrn2_)
errnum = abs(beta1)
errden = r0n2
end select
if (errden == szero) errden = sone
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
end if
beta = beta1
oldb = szero
dbar = szero
epsln = szero
phibar = beta1
rhs1 = beta1
rhs2 = szero
cs = -sone
sn = szero
itx = 0
if (do_cond) then
tnorm2 = 0.0_psb_spk_
ynorm2 = 0.0_psb_spk_
gmax = 0.0_psb_spk_
gmin = huge(sone)
end if
call psb_geaxpby(sone,r1,szero,r2,desc_a,info)
if (info == psb_success_) call psb_geaxpby(szero,x,szero,w,desc_a,info)
if (info == psb_success_) call psb_geaxpby(szero,x,szero,w1,desc_a,info)
if (info == psb_success_) call psb_geaxpby(szero,x,szero,w2,desc_a,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
if (.not.converged) then
do while (itx < itmax_)
if (abs(beta) <= epsmach) exit
itx = itx + 1
s = sone/beta
call psb_geaxpby(s,y,szero,v,desc_a,info)
if (info == psb_success_) call psb_spmm(sone,a,v,szero,y,desc_a,info,work=aux)
if (itx >= 2 .and. info == psb_success_) then
call psb_geaxpby((-beta/oldb),r1,sone,y,desc_a,info)
end if
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
alfa = psb_gedot(v,y,desc_a,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
call psb_geaxpby((-alfa/beta),r2,sone,y,desc_a,info)
if (info == psb_success_) call psb_geaxpby(sone,r2,szero,r1,desc_a,info)
if (info == psb_success_) call psb_geaxpby(sone,y,szero,r2,desc_a,info)
if (info == psb_success_) call prec%apply(r2,y,desc_a,info,work=aux)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
oldb = beta
dotv = psb_gedot(r2,y,desc_a,info)
if (dotv < -epsmach*max(sone,abs(dotv))) then
info=psb_err_pivot_too_small_
call psb_errpush(info,name,a_err='MINRES breakdown: negative (r2,z2)')
goto 9999
end if
if (dotv < szero) dotv = szero
beta = sqrt(dotv)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
if (do_cond) then
tnorm2 = tnorm2 + abs(alfa)**2 + abs(oldb)**2 + abs(beta)**2
if (itx == 1) then
gmax = abs(alfa)
gmin = gmax
end if
end if
oldeps = epsln
delta = dbar
gbar = alfa
rotx(1) = delta
roty(1) = gbar
call srot(1,rotx,1,roty,1,cs,sn)
delta = rotx(1)
gbar = -roty(1)
epsln = szero
dbar = beta
rotx(1) = epsln
roty(1) = dbar
call srot(1,rotx,1,roty,1,cs,sn)
epsln = rotx(1)
dbar = -roty(1)
gamma = gbar
betatmp = beta
call srotg(gamma,betatmp,cs,sn)
if (abs(gamma) <= epsmach) exit
phi = phibar
phibar = szero
rotx(1) = phi
roty(1) = phibar
call srot(1,rotx,1,roty,1,cs,-sn)
phi = rotx(1)
phibar = roty(1)
denom = sone/gamma
call psb_geaxpby(sone,w2,szero,w1,desc_a,info)
if (info == psb_success_) call psb_geaxpby(sone,w,szero,w2,desc_a,info)
if (info == psb_success_) call psb_geaxpby(sone,v,szero,wtmp,desc_a,info)
if (info == psb_success_) call psb_geaxpby(-oldeps,w1,sone,wtmp,desc_a,info)
if (info == psb_success_) call psb_geaxpby(-delta,w2,sone,wtmp,desc_a,info)
if (info == psb_success_) call psb_geaxpby(denom,wtmp,szero,w,desc_a,info)
if (info == psb_success_) call psb_geaxpby(phi,w,sone,x,desc_a,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
if (do_cond) then
gmax = max(gmax,abs(gamma))
gmin = min(gmin,abs(gamma))
zeta = rhs1/gamma
ynorm2 = ynorm2 + abs(zeta)**2
rhs1 = rhs2 - delta*zeta
rhs2 = -epsln*zeta
end if
select case(istop_)
case(psb_istop_ani_)
! Compute true residual only for the ANI stopping criterion.
call psb_geaxpby(sone,b,szero,res,desc_a,info)
if (info == psb_success_) call psb_spmm(-sone,a,x,sone,res,desc_a,info,work=aux)
if (info == psb_success_) rni = psb_geamax(res,desc_a,info)
if (info == psb_success_) xni = psb_geamax(x,desc_a,info)
errnum = rni
errden = ani*xni + bni
case(psb_istop_bn2_)
errnum = abs(phibar)
errden = bn2
case(psb_istop_rn2_abs_)
errnum = abs(phibar)
errden = sone
case(psb_istop_rrn2_)
errnum = abs(phibar)
errden = r0n2
end select
if (errden == szero) errden = sone
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,deps)
if (errnum <= eps*errden) then
converged = .true.
exit
end if
end do
end if
if ((.not.converged) .and. (itx >= itmax_)) then
if (debug_level >= psb_debug_ext_)&
& write(debug_unit,*) me,' ',trim(name),': MINRES reached iteration limit'
end if
if (do_cond) then
if (gmin > 0.0_psb_spk_) then
cond = gmax/gmin
else
cond = huge(sone)
end if
end if
call log_end(methdname,me,itx,itrace_,errnum,errden,deps,err=derr,iter=iter)
if (present(err)) err = derr
if (info == psb_success_) call psb_gefree(wwrk,desc_a,info)
if (info == psb_success_) deallocate(aux,stat=info)
if (info /= psb_success_) then
call psb_errpush(info,name)
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine psb_sminres_vect
linsolve/impl/psb_srgmres.f90
+133 -89
View File
@@ -97,7 +97,8 @@
! iterations
! istop - integer(optional) Input: stopping criterion, or how
! to estimate the error.
! 1: err = |r|/(|a||x|+|b|); here the iteration is
! 1: err = |r|/(|a||x|+|b|); here
! the iteration is
! stopped when |r| <= eps * (|a||x|+|b|)
! 2: err = |r|/|b|; here the iteration is
! stopped when |r| <= eps * |b|
@@ -130,7 +131,7 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
type(psb_s_vect_type) :: w, w1, xt
real(psb_spk_) :: tmp
real(psb_spk_) :: scal, gm, rti, rti1
integer(psb_ipk_) ::litmax, naux, it, k, itrace_,&
integer(psb_ipk_) ::itmax_, naux, it, k, itrace_,&
& n_row, n_col, nl
integer(psb_lpk_) :: mglob
Logical, Parameter :: exchange=.True., noexchange=.False., use_srot=.true.
@@ -172,24 +173,35 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
!
! ISTOP_ = 1: Normwise backward error, infinity norm
! ISTOP_ = 2: ||r||/||b||, 2-norm
!
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
endif
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_,psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
if (present(itmax)) then
litmax = itmax
itmax_ = itmax
else
litmax = 1000
itmax_ = 1000
endif
if (present(itrace)) then
@@ -253,12 +265,15 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
& w%get_nrows(),w1%get_nrows()
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
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)
else if (istop_ == 3) then
case(psb_istop_bn2_)
bn2 = psb_genrm2(b,desc_a,info)
case(psb_istop_rn2_abs_)
! do nothing
case(psb_istop_rrn2_)
call psb_geaxpby(sone,b,szero,v(1),desc_a,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
@@ -273,7 +288,8 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
goto 9999
end if
r0n2 = psb_genrm2(v(1),desc_a,info)
endif
end select
errnum = szero
errden = sone
deps = eps
@@ -324,27 +340,32 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
!
! check convergence
!
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
rni = psb_geamax(v(1),desc_a,info)
xni = psb_geamax(x,desc_a,info)
errnum = rni
errden = (ani*xni+bni)
else if (istop_ == 2) then
case(psb_istop_bn2_)
rni = psb_genrm2(v(1),desc_a,info)
errnum = rni
errden = bn2
else if (istop_ == 3) then
case(psb_istop_rn2_abs_)
rni = psb_genrm2(v(1),desc_a,info)
errnum = rni
errden = sone
case(psb_istop_rrn2_)
rni = psb_genrm2(v(1),desc_a,info)
errnum = rni
errden = r0n2
endif
end select
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
if ((errnum <= eps*errden).or.(itx >= litmax)) exit restart
if ((errnum <= eps*errden).or.(itx >= itmax_)) exit restart
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,deps)
@@ -360,42 +381,18 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
call prec%apply(v(i),w1,desc_a,info)
call psb_spmm(sone,a,w1,szero,w,desc_a,info,work=aux)
!
call mgs(i,h,v,w,rs,c,s,desc_a,info)
do k = 1, i
h(k,i) = psb_gedot(v(k),w,desc_a,info)
call psb_geaxpby(-h(k,i),v(k),sone,w,desc_a,info)
end do
h(i+1,i) = psb_genrm2(w,desc_a,info)
scal=sone/h(i+1,i)
call psb_geaxpby(scal,w,szero,v(i+1),desc_a,info)
do k=2,i
call srot(1,h(k-1,i),1,h(k,i),1,real(c(k-1),kind=psb_spk_),s(k-1))
enddo
rti = h(i,i)
rti1 = h(i+1,i)
call srotg(rti,rti1,tmp,s(i))
c(i) = cmplx(tmp,szero)
call srot(1,h(i,i),1,h(i+1,i),1,real(c(i),kind=psb_spk_),s(i))
h(i+1,i) = szero
call srot(1,rs(i),1,rs(i+1),1,real(c(i),kind=psb_spk_),s(i))
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
!
! build x and then compute the residual and its infinity norm
!
rst = rs
call w1%set(szero)
call strsm('l','u','n','n',i,1,sone,h,size(h,1),rst,size(rst,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rst(1:i)
do k=1, i
call psb_geaxpby(rst(k),v(k),sone,xt,desc_a,info)
end do
call prec%apply(xt,desc_a,info)
call psb_geaxpby(sone,x,sone,xt,desc_a,info)
call psb_geaxpby(sone,x,szero,xt,desc_a,info)
call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info)
call psb_geaxpby(sone,b,szero,w1,desc_a,info)
call psb_spmm(-sone,a,xt,sone,w1,desc_a,info,work=aux)
rni = psb_geamax(w1,desc_a,info)
@@ -403,8 +400,7 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
errnum = rni
errden = (ani*xni+bni)
!
else if (istop_ == 2) then
case(psb_istop_bn2_)
!
! compute the residual 2-norm as byproduct of the solution
! procedure of the least-squares problem
@@ -412,7 +408,14 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
rni = abs(rs(i+1))
errnum = rni
errden = bn2
else if (istop_ == 3) then
case(psb_istop_rn2_abs_)
rni = abs(rs(i+1))
errnum = rni
errden = sone
case(psb_istop_rrn2_)
!
! compute the residual 2-norm as byproduct of the solution
! procedure of the least-squares problem
@@ -420,28 +423,18 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
rni = abs(rs(i+1))
errnum = rni
errden = r0n2
endif
end select
if (errnum <= eps*errden) then
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
call psb_geaxpby(sone,xt,szero,x,desc_a,info)
! = x = xt
else if (istop_ == 2) then
!
! build x
!
call strsm('l','u','n','n',i,1,sone,h,size(h,1),rs,size(rs,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:i)
call w1%set(szero)
do k=1, i
call psb_geaxpby(rs(k),v(k),sone,w1,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(sone,w,sone,x,desc_a,info)
end if
case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_)
call rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info) !
end select
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,ione,errnum,errden,deps)
@@ -454,25 +447,15 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
end do inner
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
call psb_geaxpby(sone,xt,szero,x,desc_a,info)! x = xt
else if (istop_ == 2) then
!
! build x
!
call strsm('l','u','n','n',nl,1,sone,h,size(h,1),rs,size(rs,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:nl)
call w1%set(szero)
do k=1, nl
call psb_geaxpby(rs(k),v(k),sone,w1,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(sone,w,sone,x,desc_a,info)
end if
case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_)
call rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info) !
end select
if (itx >= litmax) then
if (itx >= itmax_) then
if (itrace_ > 0) then
if (mod(itx,itrace_)/=0) &
& call log_conv(methdname,me,itx,ione,errnum,errden,deps)
@@ -502,6 +485,67 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
9999 call psb_error_handler(err_act)
return
contains
!
! Advance one step the Gram-Schmidt process, and apply
! Givens rotations to keep H triangular
!
subroutine mgs(n,h,v,w,rs,c,s,desc_a,info)
real(psb_spk_) :: c(:), s(:), h(:,:), rs(:)
type(psb_s_vect_type) :: v(:), w
type(psb_desc_type) :: desc_a
real(psb_spk_) :: scal, gm, rti, rti1
real(psb_spk_) :: tmp
integer(psb_ipk_) :: info
integer(psb_ipk_) :: k,n
!
do k = 1, n
h(k,n) = psb_gedot(v(k),w,desc_a,info)
call psb_geaxpby(-h(k,n),v(k),sone,w,desc_a,info)
end do
h(n+1,n) = psb_genrm2(w,desc_a,info)
scal=sone/h(n+1,n)
call psb_geaxpby(scal,w,szero,v(n+1),desc_a,info)
do k=2,n
call srot(1,h(k-1:k-1,n),1,h(k:k,n),1,real(c(k-1),kind=psb_spk_),s(k-1))
enddo
rti = h(n,n)
rti1 = h(n+1,n)
call srotg(rti,rti1,tmp,s(n))
c(n) = cmplx(tmp,szero)
call srot(1,h(n:n,n),1,h(n+1:n+1,n),1,real(c(n),kind=psb_spk_),s(n))
call srot(1,rs(n:n),1,rs(n+1:n+1),1,real(c(n),kind=psb_spk_),s(n))
end subroutine mgs
!
! Rebuild solution X from the space V using the factor
! stored in R
!
subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info)
real(psb_spk_) :: c(:), s(:), rs(:), h(:,:)
type(psb_s_vect_type) :: v(:), w, w1, x
type(psb_desc_type) :: desc_a
class(psb_sprec_type) :: prec
integer(psb_ipk_) :: info
integer(psb_ipk_) :: k,n
!
!
! build x
!
call strsm('l','u','n','n',n,1,sone,h,size(h,1),rs,size(rs,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:n)
call w1%zero()
do k=1, n
call psb_geaxpby(rs(k),v(k),sone,w1,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(sone,w,sone,x,desc_a,info)
end subroutine rebuildx
end subroutine psb_srgmres_vect
linsolve/impl/psb_zbicg.f90
+18 -8
View File
@@ -160,19 +160,29 @@ subroutine psb_zbicg_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
endif
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if(info /= psb_success_) then
linsolve/impl/psb_zcg.F90
+22 -1
View File
@@ -159,8 +159,29 @@ subroutine psb_zcg_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_)&
linsolve/impl/psb_zcgs.f90
+22 -1
View File
@@ -154,8 +154,29 @@ Subroutine psb_zcgs_vect(a,prec,b,x,eps,desc_a,info,&
If (Present(istop)) Then
istop_ = istop
Else
istop_ = 2
istop_ = psb_get_istop_default()
Endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_) call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info)
linsolve/impl/psb_zcgstab.f90
+23 -5
View File
@@ -156,13 +156,31 @@ Subroutine psb_zcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
If (Present(istop)) Then
istop_ = istop
Else
istop_ = 2
Endif
else
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! ISTOP_ = 1: Normwise backward error, infinity norm
! ISTOP_ = 2: ||r||/||b|| norm 2
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
! = if (.not.same_type_as(x,b)) then
! = write(0,*) 'Warning: different dynamic types for X and B '
! = end if
linsolve/impl/psb_zcgstabl.f90
+22 -1
View File
@@ -172,8 +172,29 @@ Subroutine psb_zcgstabl_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
if (present(itmax)) then
itmax_ = itmax
linsolve/impl/psb_zfcg.F90
+22 -2
View File
@@ -164,9 +164,29 @@ subroutine psb_zfcg_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_)&
linsolve/impl/psb_zgcr.f90
+19 -11
View File
@@ -167,22 +167,30 @@ subroutine psb_zgcr_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
!
! ISTOP_ = 1: Normwise backward error, infinity norm
! ISTOP_ = 2: ||r||/||b||, 2-norm
!
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
endif
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_)&
& call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info)
linsolve/impl/psb_zkrylov.f90
+4 -1
View File
@@ -150,7 +150,7 @@ Subroutine psb_zkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
procedure(psb_zkryl_vect) :: psb_zbicg_vect, psb_zcgstab_vect,&
& psb_zcgs_vect
procedure(psb_zkryl_rest_vect) :: psb_zrgmres_vect, psb_zcgstabl_vect, psb_zgcr_vect
procedure(psb_zkryl_cond_vect) :: psb_zcg_vect, psb_zfcg_vect
procedure(psb_zkryl_cond_vect) :: psb_zcg_vect, psb_zfcg_vect, psb_zminres_vect
logical :: do_alloc_wrk
type(psb_ctxt_type) :: ctxt
@@ -199,6 +199,9 @@ Subroutine psb_zkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
case('RGMRES','GMRES')
call psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,irst=irst,istop=istop)
case('MINRES','PMINRES')
call psb_zminres_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,istop=istop)
case('BICGSTABL')
call psb_zcgstabl_vect(a,prec,b,x,eps,desc_a,info,&
&itmax,iter,err,itrace=itrace_,irst=irst,istop=istop)
linsolve/impl/psb_zminres.f90
+511
View File
@@ -0,0 +1,511 @@
!
! Parallel Sparse BLAS version 3.5
! (C) Copyright 2006-2018
! Salvatore Filippone
! Alfredo Buttari
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
! 3. The name of the PSBLAS group or the names of its contributors may
! not be used to endorse or promote products derived from this
! software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
!
! File: psb_zminres.f90
! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
! C C
! C References: C
! C [1] Duff, I., Marrone, M., Radicati, G., and Vittoli, C. C
! C Level 3 basic linear algebra subprograms for sparse C
! C matrices: a user level interface C
! C ACM Trans. Math. Softw., 23(3), 379-401, 1997. C
! C C
! C C
! C [2] S. Filippone, M. Colajanni C
! C PSBLAS: A library for parallel linear algebra C
! C computation on sparse matrices C
! C ACM Trans. on Math. Softw., 26(4), 527-550, Dec. 2000. C
! C C
! C [3] M. Arioli, I. Duff, M. Ruiz C
! C Stopping criteria for iterative solvers C
! C SIAM J. Matrix Anal. Appl., Vol. 13, pp. 138-144, 1992 C
! C C
! C C
! C [4] R. Barrett et al C
! C Templates for the solution of linear systems C
! C SIAM, 1993
! C C
! C C
! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
! File: psb_zminres.f90
!
! Subroutine: psb_zminres
! This subroutine implements the MINRES method.
!
!
! Arguments:
!
! a - type(psb_zspmat_type) Input: sparse matrix containing A.
! prec - class(psb_zprec_type) Input: preconditioner
! b(:) - complex Input: vector containing the
! right hand side B
! x(:) - complex Input/Output: vector containing the
! initial guess and final solution X.
! eps - real Input: Stopping tolerance; the iteration is
! stopped when the error estimate |err| <= eps
! desc_a - type(psb_desc_type). Input: The communication descriptor.
! info - integer. Output: Return code
!
! itmax - integer(optional) Input: maximum number of iterations to be
! performed.
! iter - integer(optional) Output: how many iterations have been
! performed.
! performed.
! err - real (optional) Output: error estimate on exit. If the
! denominator of the estimate is exactly
! 0, it is changed into 1.
! itrace - integer(optional) Input: print an informational message
! with the error estimate every itrace
! iterations
! istop - integer(optional) Input: stopping criterion, or how
! to estimate the error.
! 1: err = |r|/(|a||x|+|b|); here the iteration is
! stopped when |r| <= eps * (|a||x|+|b|)
! 2: err = |r|/|b|; here the iteration is
! stopped when |r| <= eps * |b|
! where r is the (preconditioned, recursive
! estimate of) residual.
!
!
subroutine psb_zminres_vect(a,prec,b,x,eps,desc_a,info,&
& itmax,iter,err,itrace,istop,cond)
use psb_base_mod
use psb_prec_mod
use psb_z_linsolve_conv_mod
use psb_linsolve_mod
implicit none
type(psb_zspmat_type), intent(in) :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(psb_zprec_type), intent(inout) :: prec
type(psb_z_vect_type), Intent(inout) :: b
type(psb_z_vect_type), Intent(inout) :: x
Real(psb_dpk_), Intent(in) :: eps
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), Optional, Intent(in) :: itmax, itrace, istop
integer(psb_ipk_), Optional, Intent(out) :: iter
Real(psb_dpk_), Optional, Intent(out) :: err, cond
! = Local data
complex(psb_dpk_), allocatable, target :: aux(:)
type(psb_z_vect_type), allocatable, target :: wwrk(:)
type(psb_z_vect_type), pointer :: y, r1, r2, v, w, w1, w2, wtmp, res
complex(psb_dpk_) :: alfa, beta, beta1, dbar, delta, denom
complex(psb_dpk_) :: epsln, gamma, gbar, oldb, oldeps, phi, phibar
complex(psb_dpk_) :: rhs1, rhs2, s, zeta, betatmp
complex(psb_dpk_) :: rotx(1), roty(1)
real(psb_dpk_) :: cs
complex(psb_dpk_) :: sn
integer(psb_ipk_) :: itmax_, istop_, naux, itx, itrace_, n_col, n_row, err_act
integer(psb_lpk_) :: mglob
integer(psb_ipk_) :: debug_level, debug_unit
type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: np, me
real(psb_dpk_) :: derr, errnum, errden, deps
logical :: do_cond
logical :: converged
real(psb_dpk_) :: epsmach
real(psb_dpk_) :: gmax, gmin, tnorm2, ynorm2
real(psb_dpk_) :: rni, xni, bni, ani, bn2, r0n2
character(len=20) :: name
character(len=*), parameter :: methdname='MINRES'
info = psb_success_
name = 'psb_zminres'
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ctxt = desc_a%get_context()
call psb_info(ctxt, me, np)
if (.not.allocated(b%v)) then
info = psb_err_invalid_vect_state_
call psb_errpush(info,name)
goto 9999
endif
if (.not.allocated(x%v)) then
info = psb_err_invalid_vect_state_
call psb_errpush(info,name)
goto 9999
endif
mglob = desc_a%get_global_rows()
n_row = desc_a%get_local_rows()
n_col = desc_a%get_local_cols()
if (present(istop)) then
istop_ = istop
else
istop_ = psb_get_istop_default()
endif
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_, psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info)
if (info == psb_success_)&
& call psb_chkvect(mglob,lone,b%get_nrows(),lone,lone,desc_a,info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='psb_chkvect on X/B')
goto 9999
end if
naux=4*n_col
allocate(aux(naux), stat=info)
if (info == psb_success_) call psb_geall(wwrk,desc_a,info,n=9_psb_ipk_)
if (info == psb_success_) call psb_geasb(wwrk,desc_a,info,mold=x%v,scratch=.true.)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
if (present(itmax)) then
itmax_ = itmax
else
itmax_ = 1000
endif
if (present(itrace)) then
itrace_ = itrace
else
itrace_ = 0
end if
y => wwrk(1)
r1 => wwrk(2)
r2 => wwrk(3)
v => wwrk(4)
w => wwrk(5)
w1 => wwrk(6)
w2 => wwrk(7)
wtmp => wwrk(8)
res => wwrk(9)
do_cond = present(cond)
epsmach = epsilon(done)
! res = b - A*x
call psb_geaxpby(zone,b,zzero,res,desc_a,info)
if (info == psb_success_) call psb_spmm(-zone,a,x,zone,res,desc_a,info,work=aux)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
select case(istop_)
case(psb_istop_ani_)
ani = psb_spnrmi(a,desc_a,info)
bni = psb_geamax(b,desc_a,info)
case(psb_istop_bn2_)
bn2 = psb_genrm2(b,desc_a,info)
case(psb_istop_rn2_abs_)
! nothing needed
case(psb_istop_rrn2_)
r0n2 = psb_genrm2(res,desc_a,info)
end select
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
errnum = zzero
errden = zone
deps = eps
converged = .false.
if ((itrace_ > 0).and.(me == 0)) call log_header(methdname)
! y = beta1 * P' * v1, with v1 the first Lanczos vector.
call psb_geaxpby(zone,res,zzero,y,desc_a,info)
if (info == psb_success_) call psb_geaxpby(zone,res,zzero,r1,desc_a,info)
if (info == psb_success_) call prec%apply(res,y,desc_a,info,work=aux)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
beta1 = sqrt(psb_gedot(res,y,desc_a,info))
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
! If beta1 is numerically zero, x is already acceptable.
if (abs(beta1) <= epsmach) then
itx = 0
converged = .true.
if (do_cond) cond = dzero
select case(istop_)
case(psb_istop_ani_)
rni = psb_geamax(res,desc_a,info)
xni = psb_geamax(x,desc_a,info)
errnum = rni
errden = ani*xni + bni
case(psb_istop_bn2_)
errnum = abs(beta1)
errden = bn2
case(psb_istop_rn2_abs_)
errnum = abs(beta1)
errden = done
case(psb_istop_rrn2_)
errnum = abs(beta1)
errden = r0n2
end select
if (errden == dzero) errden = done
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
end if
beta = beta1
oldb = zzero
dbar = zzero
epsln = zzero
phibar = beta1
rhs1 = beta1
rhs2 = zzero
cs = -done
sn = zzero
itx = 0
if (do_cond) then
tnorm2 = 0.0_psb_dpk_
ynorm2 = 0.0_psb_dpk_
gmax = 0.0_psb_dpk_
gmin = huge(done)
end if
call psb_geaxpby(zone,r1,zzero,r2,desc_a,info)
if (info == psb_success_) call psb_geaxpby(zzero,x,zzero,w,desc_a,info)
if (info == psb_success_) call psb_geaxpby(zzero,x,zzero,w1,desc_a,info)
if (info == psb_success_) call psb_geaxpby(zzero,x,zzero,w2,desc_a,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
if (.not.converged) then
do while (itx < itmax_)
if (abs(beta) <= epsmach) exit
itx = itx + 1
s = zone/beta
call psb_geaxpby(s,y,zzero,v,desc_a,info)
if (info == psb_success_) call psb_spmm(zone,a,v,zzero,y,desc_a,info,work=aux)
if (itx >= 2 .and. info == psb_success_) then
call psb_geaxpby((-beta/oldb),r1,zone,y,desc_a,info)
end if
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
alfa = psb_gedot(v,y,desc_a,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
call psb_geaxpby((-alfa/beta),r2,zone,y,desc_a,info)
if (info == psb_success_) call psb_geaxpby(zone,r2,zzero,r1,desc_a,info)
if (info == psb_success_) call psb_geaxpby(zone,y,zzero,r2,desc_a,info)
if (info == psb_success_) call prec%apply(r2,y,desc_a,info,work=aux)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
oldb = beta
beta = sqrt(psb_gedot(r2,y,desc_a,info))
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
if (do_cond) then
tnorm2 = tnorm2 + abs(alfa)**2 + abs(oldb)**2 + abs(beta)**2
if (itx == 1) then
gmax = abs(alfa)
gmin = gmax
end if
end if
oldeps = epsln
delta = dbar
gbar = alfa
rotx(1) = delta
roty(1) = gbar
call zrot(1,rotx,1,roty,1,cs,sn)
delta = rotx(1)
gbar = -roty(1)
epsln = zzero
dbar = beta
rotx(1) = epsln
roty(1) = dbar
call zrot(1,rotx,1,roty,1,cs,sn)
epsln = rotx(1)
dbar = -roty(1)
gamma = gbar
betatmp = beta
call zrotg(gamma,betatmp,cs,sn)
if (abs(gamma) <= epsmach) exit
phi = phibar
phibar = zzero
rotx(1) = phi
roty(1) = phibar
call zrot(1,rotx,1,roty,1,cs,-sn)
phi = rotx(1)
phibar = roty(1)
denom = zone/gamma
call psb_geaxpby(zone,w2,zzero,w1,desc_a,info)
if (info == psb_success_) call psb_geaxpby(zone,w,zzero,w2,desc_a,info)
if (info == psb_success_) call psb_geaxpby(zone,v,zzero,wtmp,desc_a,info)
if (info == psb_success_) call psb_geaxpby(-oldeps,w1,zone,wtmp,desc_a,info)
if (info == psb_success_) call psb_geaxpby(-delta,w2,zone,wtmp,desc_a,info)
if (info == psb_success_) call psb_geaxpby(denom,wtmp,zzero,w,desc_a,info)
if (info == psb_success_) call psb_geaxpby(phi,w,zone,x,desc_a,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
if (do_cond) then
gmax = max(gmax,abs(gamma))
gmin = min(gmin,abs(gamma))
zeta = rhs1/gamma
ynorm2 = ynorm2 + abs(zeta)**2
rhs1 = rhs2 - delta*zeta
rhs2 = -epsln*zeta
end if
select case(istop_)
case(psb_istop_ani_)
! Compute true residual only for the ANI stopping criterion.
call psb_geaxpby(zone,b,zzero,res,desc_a,info)
if (info == psb_success_) call psb_spmm(-zone,a,x,zone,res,desc_a,info,work=aux)
if (info == psb_success_) rni = psb_geamax(res,desc_a,info)
if (info == psb_success_) xni = psb_geamax(x,desc_a,info)
errnum = rni
errden = ani*xni + bni
case(psb_istop_bn2_)
errnum = abs(phibar)
errden = bn2
case(psb_istop_rn2_abs_)
errnum = abs(phibar)
errden = done
case(psb_istop_rrn2_)
errnum = abs(phibar)
errden = r0n2
end select
if (errden == dzero) errden = done
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,deps)
if (errnum <= eps*errden) then
converged = .true.
exit
end if
end do
end if
if ((.not.converged) .and. (itx >= itmax_)) then
if (debug_level >= psb_debug_ext_)&
& write(debug_unit,*) me,' ',trim(name),': MINRES reached iteration limit'
end if
if (do_cond) then
if (gmin > 0.0_psb_dpk_) then
cond = gmax/gmin
else
cond = huge(done)
end if
end if
call log_end(methdname,me,itx,itrace_,errnum,errden,deps,err=derr,iter=iter)
if (present(err)) err = derr
if (info == psb_success_) call psb_gefree(wwrk,desc_a,info)
if (info == psb_success_) deallocate(aux,stat=info)
if (info /= psb_success_) then
call psb_errpush(info,name)
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine psb_zminres_vect
linsolve/impl/psb_zrgmres.f90
+133 -89
View File
@@ -97,7 +97,8 @@
! iterations
! istop - integer(optional) Input: stopping criterion, or how
! to estimate the error.
! 1: err = |r|/(|a||x|+|b|); here the iteration is
! 1: err = |r|/(|a||x|+|b|); here
! the iteration is
! stopped when |r| <= eps * (|a||x|+|b|)
! 2: err = |r|/|b|; here the iteration is
! stopped when |r| <= eps * |b|
@@ -130,7 +131,7 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
type(psb_z_vect_type) :: w, w1, xt
real(psb_dpk_) :: tmp
complex(psb_dpk_) :: scal, gm, rti, rti1
integer(psb_ipk_) ::litmax, naux, it, k, itrace_,&
integer(psb_ipk_) ::itmax_, naux, it, k, itrace_,&
& n_row, n_col, nl
integer(psb_lpk_) :: mglob
Logical, Parameter :: exchange=.True., noexchange=.False., use_srot=.true.
@@ -172,24 +173,35 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
if (present(istop)) then
istop_ = istop
else
istop_ = 2
istop_ = psb_get_istop_default()
endif
!
! ISTOP_ = 1: Normwise backward error, infinity norm
! ISTOP_ = 2: ||r||/||b||, 2-norm
!
if ((istop_ < 1 ).or.(istop_ > 2 ) ) then
if (.not.psb_is_valid_istop(istop_)) then
info=psb_err_invalid_istop_
err=info
call psb_errpush(info,name,i_err=(/istop_/))
goto 9999
endif
end if
!
! istop_ = 1: normwise backward error, infinity norm
! istop_ = 2: ||r||/||b|| norm 2
!
select case(istop_)
case(psb_istop_ani_,psb_istop_bn2_,&
& psb_istop_rn2_abs_,psb_istop_rrn2_)
! nothing needed
case default
! should never get here
info=psb_err_internal_error_
err=info
call psb_errpush(info,name,a_err="invalid istop_")
goto 9999
end select
if (present(itmax)) then
litmax = itmax
itmax_ = itmax
else
litmax = 1000
itmax_ = 1000
endif
if (present(itrace)) then
@@ -253,12 +265,15 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
& w%get_nrows(),w1%get_nrows()
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
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)
else if (istop_ == 3) then
case(psb_istop_bn2_)
bn2 = psb_genrm2(b,desc_a,info)
case(psb_istop_rn2_abs_)
! do nothing
case(psb_istop_rrn2_)
call psb_geaxpby(zone,b,zzero,v(1),desc_a,info)
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
@@ -273,7 +288,8 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
goto 9999
end if
r0n2 = psb_genrm2(v(1),desc_a,info)
endif
end select
errnum = zzero
errden = zone
deps = eps
@@ -324,27 +340,32 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
!
! check convergence
!
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
rni = psb_geamax(v(1),desc_a,info)
xni = psb_geamax(x,desc_a,info)
errnum = rni
errden = (ani*xni+bni)
else if (istop_ == 2) then
case(psb_istop_bn2_)
rni = psb_genrm2(v(1),desc_a,info)
errnum = rni
errden = bn2
else if (istop_ == 3) then
case(psb_istop_rn2_abs_)
rni = psb_genrm2(v(1),desc_a,info)
errnum = rni
errden = done
case(psb_istop_rrn2_)
rni = psb_genrm2(v(1),desc_a,info)
errnum = rni
errden = r0n2
endif
end select
if (info /= psb_success_) then
info=psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
if ((errnum <= eps*errden).or.(itx >= litmax)) exit restart
if ((errnum <= eps*errden).or.(itx >= itmax_)) exit restart
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,itrace_,errnum,errden,deps)
@@ -360,42 +381,18 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
call prec%apply(v(i),w1,desc_a,info)
call psb_spmm(zone,a,w1,zzero,w,desc_a,info,work=aux)
!
call mgs(i,h,v,w,rs,c,s,desc_a,info)
do k = 1, i
h(k,i) = psb_gedot(v(k),w,desc_a,info)
call psb_geaxpby(-h(k,i),v(k),zone,w,desc_a,info)
end do
h(i+1,i) = psb_genrm2(w,desc_a,info)
scal=zone/h(i+1,i)
call psb_geaxpby(scal,w,zzero,v(i+1),desc_a,info)
do k=2,i
call zrot(1,h(k-1,i),1,h(k,i),1,real(c(k-1),kind=psb_dpk_),s(k-1))
enddo
rti = h(i,i)
rti1 = h(i+1,i)
call zrotg(rti,rti1,tmp,s(i))
c(i) = cmplx(tmp,dzero)
call zrot(1,h(i,i),1,h(i+1,i),1,real(c(i),kind=psb_dpk_),s(i))
h(i+1,i) = zzero
call zrot(1,rs(i),1,rs(i+1),1,real(c(i),kind=psb_dpk_),s(i))
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
!
! build x and then compute the residual and its infinity norm
!
rst = rs
call w1%set(zzero)
call ztrsm('l','u','n','n',i,1,zone,h,size(h,1),rst,size(rst,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rst(1:i)
do k=1, i
call psb_geaxpby(rst(k),v(k),zone,xt,desc_a,info)
end do
call prec%apply(xt,desc_a,info)
call psb_geaxpby(zone,x,zone,xt,desc_a,info)
call psb_geaxpby(zone,x,zzero,xt,desc_a,info)
call rebuildx(i,h,v,w,w1,xt,rst,c,s,prec,desc_a,info)
call psb_geaxpby(zone,b,zzero,w1,desc_a,info)
call psb_spmm(-zone,a,xt,zone,w1,desc_a,info,work=aux)
rni = psb_geamax(w1,desc_a,info)
@@ -403,8 +400,7 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
errnum = rni
errden = (ani*xni+bni)
!
else if (istop_ == 2) then
case(psb_istop_bn2_)
!
! compute the residual 2-norm as byproduct of the solution
! procedure of the least-squares problem
@@ -412,7 +408,14 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
rni = abs(rs(i+1))
errnum = rni
errden = bn2
else if (istop_ == 3) then
case(psb_istop_rn2_abs_)
rni = abs(rs(i+1))
errnum = rni
errden = done
case(psb_istop_rrn2_)
!
! compute the residual 2-norm as byproduct of the solution
! procedure of the least-squares problem
@@ -420,28 +423,18 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
rni = abs(rs(i+1))
errnum = rni
errden = r0n2
endif
end select
if (errnum <= eps*errden) then
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
call psb_geaxpby(zone,xt,zzero,x,desc_a,info)
! = x = xt
else if (istop_ == 2) then
!
! build x
!
call ztrsm('l','u','n','n',i,1,zone,h,size(h,1),rs,size(rs,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:i)
call w1%set(zzero)
do k=1, i
call psb_geaxpby(rs(k),v(k),zone,w1,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(zone,w,zone,x,desc_a,info)
end if
case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_)
call rebuildx(i,h,v,w,w1,x,rs,c,s,prec,desc_a,info) !
end select
if (itrace_ > 0) &
& call log_conv(methdname,me,itx,ione,errnum,errden,deps)
@@ -454,25 +447,15 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
end do inner
if (istop_ == 1) then
select case(istop_)
case(psb_istop_ani_)
call psb_geaxpby(zone,xt,zzero,x,desc_a,info)! x = xt
else if (istop_ == 2) then
!
! build x
!
call ztrsm('l','u','n','n',nl,1,zone,h,size(h,1),rs,size(rs,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:nl)
call w1%set(zzero)
do k=1, nl
call psb_geaxpby(rs(k),v(k),zone,w1,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(zone,w,zone,x,desc_a,info)
end if
case(psb_istop_bn2_, psb_istop_rn2_abs_,psb_istop_rrn2_)
call rebuildx(nl,h,v,w,w1,x,rs,c,s,prec,desc_a,info) !
end select
if (itx >= litmax) then
if (itx >= itmax_) then
if (itrace_ > 0) then
if (mod(itx,itrace_)/=0) &
& call log_conv(methdname,me,itx,ione,errnum,errden,deps)
@@ -502,6 +485,67 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
9999 call psb_error_handler(err_act)
return
contains
!
! Advance one step the Gram-Schmidt process, and apply
! Givens rotations to keep H triangular
!
subroutine mgs(n,h,v,w,rs,c,s,desc_a,info)
complex(psb_dpk_) :: c(:), s(:), h(:,:), rs(:)
type(psb_z_vect_type) :: v(:), w
type(psb_desc_type) :: desc_a
complex(psb_dpk_) :: scal, gm, rti, rti1
real(psb_dpk_) :: tmp
integer(psb_ipk_) :: info
integer(psb_ipk_) :: k,n
!
do k = 1, n
h(k,n) = psb_gedot(v(k),w,desc_a,info)
call psb_geaxpby(-h(k,n),v(k),zone,w,desc_a,info)
end do
h(n+1,n) = psb_genrm2(w,desc_a,info)
scal=zone/h(n+1,n)
call psb_geaxpby(scal,w,zzero,v(n+1),desc_a,info)
do k=2,n
call zrot(1,h(k-1:k-1,n),1,h(k:k,n),1,real(c(k-1),kind=psb_dpk_),s(k-1))
enddo
rti = h(n,n)
rti1 = h(n+1,n)
call zrotg(rti,rti1,tmp,s(n))
c(n) = cmplx(tmp,dzero)
call zrot(1,h(n:n,n),1,h(n+1:n+1,n),1,real(c(n),kind=psb_dpk_),s(n))
call zrot(1,rs(n:n),1,rs(n+1:n+1),1,real(c(n),kind=psb_dpk_),s(n))
end subroutine mgs
!
! Rebuild solution X from the space V using the factor
! stored in R
!
subroutine rebuildx(n,h,v,w,w1,x,rs,c,s,prec,desc_a,info)
complex(psb_dpk_) :: c(:), s(:), rs(:), h(:,:)
type(psb_z_vect_type) :: v(:), w, w1, x
type(psb_desc_type) :: desc_a
class(psb_zprec_type) :: prec
integer(psb_ipk_) :: info
integer(psb_ipk_) :: k,n
!
!
! build x
!
call ztrsm('l','u','n','n',n,1,zone,h,size(h,1),rs,size(rs,1))
if (debug_level >= psb_debug_ext_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Rebuild x-> RS:',rs(1:n)
call w1%zero()
do k=1, n
call psb_geaxpby(rs(k),v(k),zone,w1,desc_a,info)
end do
call prec%apply(w1,w,desc_a,info)
call psb_geaxpby(zone,w,zone,x,desc_a,info)
end subroutine rebuildx
end subroutine psb_zrgmres_vect
linsolve/psb_base_linsolve_conv_mod.f90
+32 -1
View File
@@ -41,10 +41,20 @@ Module psb_base_linsolve_conv_mod
module procedure psb_d_end_conv
end interface
!
integer(psb_ipk_), parameter :: psb_istop_min_ = 1
integer(psb_ipk_), parameter :: psb_istop_ani_ = 1
integer(psb_ipk_), parameter :: psb_istop_bn2_ = 2
integer(psb_ipk_), parameter :: psb_istop_rn2_abs_ = 3
integer(psb_ipk_), parameter :: psb_istop_rrn2_ = 4
integer(psb_ipk_), parameter :: psb_istop_scbn2_ = 5
integer(psb_ipk_), parameter :: psb_istop_max_ = 3
! Fields in controls and values
integer(psb_ipk_), parameter :: psb_ik_bni_=1, psb_ik_rni_=2, psb_ik_ani_=3
integer(psb_ipk_), parameter :: psb_ik_xni_=4, psb_ik_bn2_=5, psb_ik_r0n2_=6
integer(psb_ipk_), parameter :: psb_ik_xn2_=7, psb_ik_errnum_=8, psb_ik_errden_=9
integer(psb_ipk_), parameter :: psb_ik_eps_=10, psb_ik_rn2_=11
integer(psb_ipk_), parameter :: psb_ik_eps_=10, psb_ik_rn2_=11, psb_ik_rn2_abs_=12
integer(psb_ipk_), parameter :: psb_ik_stopc_=1, psb_ik_trace_=2, psb_ik_itmax_=3
integer(psb_ipk_), parameter :: psb_ik_ivsz_=16
type psb_itconv_type
@@ -52,8 +62,29 @@ Module psb_base_linsolve_conv_mod
real(psb_dpk_) :: values(psb_ik_ivsz_)
end type psb_itconv_type
integer(psb_ipk_), save :: psb_istop_default = psb_istop_bn2_
contains
function psb_is_valid_istop(istop) result(res)
integer(psb_ipk_) :: istop
logical :: res
res = ((psb_istop_min_<=istop).and.(istop<=psb_istop_max_))
end function psb_is_valid_istop
function psb_get_istop_default() result(res)
integer(psb_ipk_) :: res
res = psb_istop_default
end function psb_get_istop_default
subroutine psb_set_istop_default(val)
integer(psb_ipk_) :: val
if ((psb_istop_min_<=val).and.(val<=psb_istop_max_)) &
& psb_istop_default = val
end subroutine psb_set_istop_default
subroutine log_header(methdname)
!use psb_base_mod
implicit none
test/pdegen/psb_d_pde2d.F90
+3 -3
View File
@@ -467,8 +467,8 @@ contains
! compute gridpoint coordinates
call idx2ijk(ix,iy,glob_row,idim,idim)
! x, y coordinates
x = (ix-1)*deltah
y = (iy-1)*deltah
x = (ix)*deltah
y = (iy)*deltah
zt(k) = f_(x,y)
! internal point: build discretization
@@ -778,7 +778,7 @@ program psb_d_pde2d
& desc_a,info,itmax=itmax,iter=iter,&
& err=err,itrace=itrace,&
& istop=istopc)
case('BICGSTAB','BICGSTABL','BICG','CG','CGS','FCG','GCR','RGMRES')
case('BICGSTAB','BICGSTABL','BICG','CG','CGS','FCG','GCR','RGMRES','MINRES')
call psb_krylov(kmethd,a,prec,bv,xxv,eps,&
& desc_a,info,itmax=itmax,iter=iter,err=err,itrace=itrace,&
& istop=istopc,irst=irst)
test/pdegen/psb_d_pde3d.F90
+4 -4
View File
@@ -496,9 +496,9 @@ contains
! compute gridpoint coordinates
call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim)
! x, y, z coordinates
x = (ix-1)*deltah
y = (iy-1)*deltah
z = (iz-1)*deltah
x = (ix)*deltah
y = (iy)*deltah
z = (iz)*deltah
zt(k) = f_(x,y,z)
! internal point: build discretization
!
@@ -834,7 +834,7 @@ program psb_d_pde3d
& desc_a,info,itmax=itmax,iter=iter,&
& err=err,itrace=itrace,&
& istop=istopc)
case('BICGSTAB','BICGSTABL','BICG','CG','CGS','FCG','GCR','RGMRES')
case('BICGSTAB','BICGSTABL','BICG','CG','CGS','FCG','GCR','RGMRES','MINRES')
call psb_krylov(kmethd,a,prec,bv,xxv,eps,&
& desc_a,info,itmax=itmax,iter=iter,err=err,itrace=itrace,&
& istop=istopc,irst=irst)
test/pdegen/psb_s_pde2d.F90
+3 -3
View File
@@ -467,8 +467,8 @@ contains
! compute gridpoint coordinates
call idx2ijk(ix,iy,glob_row,idim,idim)
! x, y coordinates
x = (ix-1)*deltah
y = (iy-1)*deltah
x = (ix)*deltah
y = (iy)*deltah
zt(k) = f_(x,y)
! internal point: build discretization
@@ -778,7 +778,7 @@ program psb_s_pde2d
& desc_a,info,itmax=itmax,iter=iter,&
& err=err,itrace=itrace,&
& istop=istopc)
case('BICGSTAB','BICGSTABL','BICG','CG','CGS','FCG','GCR','RGMRES')
case('BICGSTAB','BICGSTABL','BICG','CG','CGS','FCG','GCR','RGMRES','MINRES')
call psb_krylov(kmethd,a,prec,bv,xxv,eps,&
& desc_a,info,itmax=itmax,iter=iter,err=err,itrace=itrace,&
& istop=istopc,irst=irst)
test/pdegen/psb_s_pde3d.F90
+4 -4
View File
@@ -496,9 +496,9 @@ contains
! compute gridpoint coordinates
call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim)
! x, y, z coordinates
x = (ix-1)*deltah
y = (iy-1)*deltah
z = (iz-1)*deltah
x = (ix)*deltah
y = (iy)*deltah
z = (iz)*deltah
zt(k) = f_(x,y,z)
! internal point: build discretization
!
@@ -834,7 +834,7 @@ program psb_s_pde3d
& desc_a,info,itmax=itmax,iter=iter,&
& err=err,itrace=itrace,&
& istop=istopc)
case('BICGSTAB','BICGSTABL','BICG','CG','CGS','FCG','GCR','RGMRES')
case('BICGSTAB','BICGSTABL','BICG','CG','CGS','FCG','GCR','RGMRES','MINRES')
call psb_krylov(kmethd,a,prec,bv,xxv,eps,&
& desc_a,info,itmax=itmax,iter=iter,err=err,itrace=itrace,&
& istop=istopc,irst=irst)

Some files were not shown because too many files have changed in this diff Show More