diff --git a/base/modules/Makefile b/base/modules/Makefile index 353f64211..c43ab28e1 100644 --- a/base/modules/Makefile +++ b/base/modules/Makefile @@ -196,7 +196,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 comm/comm_schemes/psb_comm_schemes_mod.o: psb_const_mod.o @@ -463,4 +463,4 @@ clean: /bin/rm -f $(MODULES) $(OBJS) $(MPFOBJS) *$(.mod) veryclean: clean - /bin/rm -f *.h \ No newline at end of file + /bin/rm -f *.h diff --git a/base/modules/desc/psb_desc_mod.F90 b/base/modules/desc/psb_desc_mod.F90 index 4b109465f..d48ef4818 100644 --- a/base/modules/desc/psb_desc_mod.F90 +++ b/base/modules/desc/psb_desc_mod.F90 @@ -685,7 +685,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 diff --git a/config/pac.m4 b/config/pac.m4 index 9cb556476..1b2d7b757 100644 --- a/config/pac.m4 +++ b/config/pac.m4 @@ -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 ] ) diff --git a/configure b/configure index e8a214af5..d81424b36 100755 --- a/configure +++ b/configure @@ -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} diff --git a/configure.ac b/configure.ac index 1116417a9..338d05be8 100644 --- a/configure.ac +++ b/configure.ac @@ -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} diff --git a/cuda/Makefile b/cuda/Makefile index ee4757ce0..61ba47271 100644 --- a/cuda/Makefile +++ b/cuda/Makefile @@ -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 \ diff --git a/cuda/cuda_util.c b/cuda/cuda_util.c index 707c88de1..cacdb1b37 100644 --- a/cuda/cuda_util.c +++ b/cuda/cuda_util.c @@ -338,6 +338,7 @@ int getGPUMemoryClockRate() #endif return(count); } + int getGPUWarpSize() { int count=0; if (prop!=NULL) diff --git a/ext/impl/psb_c_hll_csmv.f90 b/ext/impl/psb_c_hll_csmv.f90 index ce597e2e3..f50104826 100644 --- a/ext/impl/psb_c_hll_csmv.f90 +++ b/ext/impl/psb_c_hll_csmv.f90 @@ -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_ diff --git a/ext/impl/psb_d_hll_csmv.f90 b/ext/impl/psb_d_hll_csmv.f90 index 52816fd78..f3ccb49b7 100644 --- a/ext/impl/psb_d_hll_csmv.f90 +++ b/ext/impl/psb_d_hll_csmv.f90 @@ -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_ diff --git a/ext/impl/psb_s_hll_csmv.f90 b/ext/impl/psb_s_hll_csmv.f90 index a3ead4c5c..8c8e4a288 100644 --- a/ext/impl/psb_s_hll_csmv.f90 +++ b/ext/impl/psb_s_hll_csmv.f90 @@ -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_ diff --git a/ext/impl/psb_z_hll_csmv.f90 b/ext/impl/psb_z_hll_csmv.f90 index cf871d601..808395f54 100644 --- a/ext/impl/psb_z_hll_csmv.f90 +++ b/ext/impl/psb_z_hll_csmv.f90 @@ -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_ diff --git a/linsolve/impl/CMakeLists.txt b/linsolve/impl/CMakeLists.txt index 9f2fe2d5d..d9760bb2f 100644 --- a/linsolve/impl/CMakeLists.txt +++ b/linsolve/impl/CMakeLists.txt @@ -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}) diff --git a/linsolve/impl/Makefile b/linsolve/impl/Makefile index 62444af7f..66ff22433 100644 --- a/linsolve/impl/Makefile +++ b/linsolve/impl/Makefile @@ -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= diff --git a/linsolve/impl/psb_cbicg.f90 b/linsolve/impl/psb_cbicg.f90 index 15a1aa0a2..e346b3ff0 100644 --- a/linsolve/impl/psb_cbicg.f90 +++ b/linsolve/impl/psb_cbicg.f90 @@ -159,19 +159,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 + 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_ = 2: ||r||/||b|| norm 2 ! - - if ((istop_ < 1 ).or.(istop_ > 2 ) ) then - info=psb_err_invalid_istop_ + 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,i_err=(/istop_/)) + call psb_errpush(info,name,a_err="invalid istop_") goto 9999 - endif + end select call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info) if(info /= psb_success_) then diff --git a/linsolve/impl/psb_ccg.F90 b/linsolve/impl/psb_ccg.F90 index c4cca5937..b703bb1c1 100644 --- a/linsolve/impl/psb_ccg.F90 +++ b/linsolve/impl/psb_ccg.F90 @@ -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_)& diff --git a/linsolve/impl/psb_ccgs.f90 b/linsolve/impl/psb_ccgs.f90 index 2db1d65ce..a29a74e99 100644 --- a/linsolve/impl/psb_ccgs.f90 +++ b/linsolve/impl/psb_ccgs.f90 @@ -153,8 +153,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) diff --git a/linsolve/impl/psb_ccgstab.f90 b/linsolve/impl/psb_ccgstab.f90 index 1688024a8..b4bb17d22 100644 --- a/linsolve/impl/psb_ccgstab.f90 +++ b/linsolve/impl/psb_ccgstab.f90 @@ -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 diff --git a/linsolve/impl/psb_ccgstabl.f90 b/linsolve/impl/psb_ccgstabl.f90 index ea66523ed..b11a7910b 100644 --- a/linsolve/impl/psb_ccgstabl.f90 +++ b/linsolve/impl/psb_ccgstabl.f90 @@ -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 diff --git a/linsolve/impl/psb_cfcg.F90 b/linsolve/impl/psb_cfcg.F90 index 910b611de..3c0b02591 100644 --- a/linsolve/impl/psb_cfcg.F90 +++ b/linsolve/impl/psb_cfcg.F90 @@ -163,9 +163,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_)& diff --git a/linsolve/impl/psb_cgcr.f90 b/linsolve/impl/psb_cgcr.f90 index 99865aec4..f553c5642 100644 --- a/linsolve/impl/psb_cgcr.f90 +++ b/linsolve/impl/psb_cgcr.f90 @@ -166,22 +166,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) diff --git a/linsolve/impl/psb_ckrylov.f90 b/linsolve/impl/psb_ckrylov.f90 index 22712da05..308e62aea 100644 --- a/linsolve/impl/psb_ckrylov.f90 +++ b/linsolve/impl/psb_ckrylov.f90 @@ -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) diff --git a/linsolve/impl/psb_cminres.f90 b/linsolve/impl/psb_cminres.f90 new file mode 100644 index 000000000..f4b1f8c77 --- /dev/null +++ b/linsolve/impl/psb_cminres.f90 @@ -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 diff --git a/linsolve/impl/psb_crgmres.f90 b/linsolve/impl/psb_crgmres.f90 index 214a89b39..729fc0528 100644 --- a/linsolve/impl/psb_crgmres.f90 +++ b/linsolve/impl/psb_crgmres.f90 @@ -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| @@ -129,7 +130,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, 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. @@ -171,24 +172,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 @@ -247,12 +259,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_ @@ -267,7 +282,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 @@ -318,27 +334,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) @@ -354,42 +375,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) ! + 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) rni = psb_geamax(w1,desc_a,info) @@ -397,8 +394,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 @@ -406,7 +402,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 @@ -414,28 +417,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) @@ -448,25 +441,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) @@ -496,6 +479,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 diff --git a/linsolve/impl/psb_dbicg.f90 b/linsolve/impl/psb_dbicg.f90 index 1ac471995..dad868b43 100644 --- a/linsolve/impl/psb_dbicg.f90 +++ b/linsolve/impl/psb_dbicg.f90 @@ -159,19 +159,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 + 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_ = 2: ||r||/||b|| norm 2 ! - - if ((istop_ < 1 ).or.(istop_ > 2 ) ) then - info=psb_err_invalid_istop_ + 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,i_err=(/istop_/)) + call psb_errpush(info,name,a_err="invalid istop_") goto 9999 - endif + end select call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info) if(info /= psb_success_) then diff --git a/linsolve/impl/psb_dcg.F90 b/linsolve/impl/psb_dcg.F90 index e7f9cf941..fecdeaece 100644 --- a/linsolve/impl/psb_dcg.F90 +++ b/linsolve/impl/psb_dcg.F90 @@ -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_)& diff --git a/linsolve/impl/psb_dcgs.f90 b/linsolve/impl/psb_dcgs.f90 index ad4849b8e..de6fc164d 100644 --- a/linsolve/impl/psb_dcgs.f90 +++ b/linsolve/impl/psb_dcgs.f90 @@ -153,8 +153,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) diff --git a/linsolve/impl/psb_dcgstab.f90 b/linsolve/impl/psb_dcgstab.f90 index 4b9e426ff..1cc395a66 100644 --- a/linsolve/impl/psb_dcgstab.f90 +++ b/linsolve/impl/psb_dcgstab.f90 @@ -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 diff --git a/linsolve/impl/psb_dcgstabl.f90 b/linsolve/impl/psb_dcgstabl.f90 index 750097833..f09529f0f 100644 --- a/linsolve/impl/psb_dcgstabl.f90 +++ b/linsolve/impl/psb_dcgstabl.f90 @@ -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 diff --git a/linsolve/impl/psb_dfcg.F90 b/linsolve/impl/psb_dfcg.F90 index 3b7dcea4f..e5562b4ed 100644 --- a/linsolve/impl/psb_dfcg.F90 +++ b/linsolve/impl/psb_dfcg.F90 @@ -163,9 +163,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_)& diff --git a/linsolve/impl/psb_dgcr.f90 b/linsolve/impl/psb_dgcr.f90 index c9bc8e6d1..941413377 100644 --- a/linsolve/impl/psb_dgcr.f90 +++ b/linsolve/impl/psb_dgcr.f90 @@ -166,22 +166,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) diff --git a/linsolve/impl/psb_dkrylov.f90 b/linsolve/impl/psb_dkrylov.f90 index 5963f5a8a..37daca4f9 100644 --- a/linsolve/impl/psb_dkrylov.f90 +++ b/linsolve/impl/psb_dkrylov.f90 @@ -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) diff --git a/linsolve/impl/psb_dminres.f90 b/linsolve/impl/psb_dminres.f90 new file mode 100644 index 000000000..a3e1675a1 --- /dev/null +++ b/linsolve/impl/psb_dminres.f90 @@ -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 diff --git a/linsolve/impl/psb_drgmres.f90 b/linsolve/impl/psb_drgmres.f90 index 584b93876..b86af5922 100644 --- a/linsolve/impl/psb_drgmres.f90 +++ b/linsolve/impl/psb_drgmres.f90 @@ -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| @@ -129,7 +130,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, 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. @@ -171,24 +172,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 @@ -247,12 +259,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_ @@ -267,7 +282,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 @@ -318,27 +334,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) @@ -354,42 +375,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) ! + 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) rni = psb_geamax(w1,desc_a,info) @@ -397,8 +394,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 @@ -406,7 +402,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 @@ -414,28 +417,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) @@ -448,25 +441,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) @@ -496,6 +479,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 diff --git a/linsolve/impl/psb_sbicg.f90 b/linsolve/impl/psb_sbicg.f90 index bb687f037..7cda3ba95 100644 --- a/linsolve/impl/psb_sbicg.f90 +++ b/linsolve/impl/psb_sbicg.f90 @@ -159,19 +159,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 + 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_ = 2: ||r||/||b|| norm 2 ! - - if ((istop_ < 1 ).or.(istop_ > 2 ) ) then - info=psb_err_invalid_istop_ + 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,i_err=(/istop_/)) + call psb_errpush(info,name,a_err="invalid istop_") goto 9999 - endif + end select call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info) if(info /= psb_success_) then diff --git a/linsolve/impl/psb_scg.F90 b/linsolve/impl/psb_scg.F90 index 198537eb1..037d42e6b 100644 --- a/linsolve/impl/psb_scg.F90 +++ b/linsolve/impl/psb_scg.F90 @@ -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_)& diff --git a/linsolve/impl/psb_scgs.f90 b/linsolve/impl/psb_scgs.f90 index f6eb345e6..74d314252 100644 --- a/linsolve/impl/psb_scgs.f90 +++ b/linsolve/impl/psb_scgs.f90 @@ -153,8 +153,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) diff --git a/linsolve/impl/psb_scgstab.f90 b/linsolve/impl/psb_scgstab.f90 index 2b2bdcfa3..c97e0ae12 100644 --- a/linsolve/impl/psb_scgstab.f90 +++ b/linsolve/impl/psb_scgstab.f90 @@ -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 diff --git a/linsolve/impl/psb_scgstabl.f90 b/linsolve/impl/psb_scgstabl.f90 index 125e3b92d..4dcf58358 100644 --- a/linsolve/impl/psb_scgstabl.f90 +++ b/linsolve/impl/psb_scgstabl.f90 @@ -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 diff --git a/linsolve/impl/psb_sfcg.F90 b/linsolve/impl/psb_sfcg.F90 index 940aa21ef..1c071cc23 100644 --- a/linsolve/impl/psb_sfcg.F90 +++ b/linsolve/impl/psb_sfcg.F90 @@ -163,9 +163,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_)& diff --git a/linsolve/impl/psb_sgcr.f90 b/linsolve/impl/psb_sgcr.f90 index 6c45d44be..40429d366 100644 --- a/linsolve/impl/psb_sgcr.f90 +++ b/linsolve/impl/psb_sgcr.f90 @@ -166,22 +166,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) diff --git a/linsolve/impl/psb_skrylov.f90 b/linsolve/impl/psb_skrylov.f90 index e24d7d05f..cce35cf9a 100644 --- a/linsolve/impl/psb_skrylov.f90 +++ b/linsolve/impl/psb_skrylov.f90 @@ -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) diff --git a/linsolve/impl/psb_sminres.f90 b/linsolve/impl/psb_sminres.f90 new file mode 100644 index 000000000..7696e323a --- /dev/null +++ b/linsolve/impl/psb_sminres.f90 @@ -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 diff --git a/linsolve/impl/psb_srgmres.f90 b/linsolve/impl/psb_srgmres.f90 index a0db5f698..320d3166d 100644 --- a/linsolve/impl/psb_srgmres.f90 +++ b/linsolve/impl/psb_srgmres.f90 @@ -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| @@ -129,7 +130,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, 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. @@ -171,24 +172,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 @@ -247,12 +259,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_ @@ -267,7 +282,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 @@ -318,27 +334,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) @@ -354,42 +375,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) ! + 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) rni = psb_geamax(w1,desc_a,info) @@ -397,8 +394,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 @@ -406,7 +402,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 @@ -414,28 +417,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) @@ -448,25 +441,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) @@ -496,6 +479,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 diff --git a/linsolve/impl/psb_zbicg.f90 b/linsolve/impl/psb_zbicg.f90 index 19ba88001..017ee1073 100644 --- a/linsolve/impl/psb_zbicg.f90 +++ b/linsolve/impl/psb_zbicg.f90 @@ -159,19 +159,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 + 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_ = 2: ||r||/||b|| norm 2 ! - - if ((istop_ < 1 ).or.(istop_ > 2 ) ) then - info=psb_err_invalid_istop_ + 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,i_err=(/istop_/)) + call psb_errpush(info,name,a_err="invalid istop_") goto 9999 - endif + end select call psb_chkvect(mglob,lone,x%get_nrows(),lone,lone,desc_a,info) if(info /= psb_success_) then diff --git a/linsolve/impl/psb_zcg.F90 b/linsolve/impl/psb_zcg.F90 index dbc6ea2ee..9835f0c1a 100644 --- a/linsolve/impl/psb_zcg.F90 +++ b/linsolve/impl/psb_zcg.F90 @@ -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_)& diff --git a/linsolve/impl/psb_zcgs.f90 b/linsolve/impl/psb_zcgs.f90 index 88d0d927d..a67ecebcf 100644 --- a/linsolve/impl/psb_zcgs.f90 +++ b/linsolve/impl/psb_zcgs.f90 @@ -153,8 +153,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) diff --git a/linsolve/impl/psb_zcgstab.f90 b/linsolve/impl/psb_zcgstab.f90 index 16bb5ddc5..e143d7f7a 100644 --- a/linsolve/impl/psb_zcgstab.f90 +++ b/linsolve/impl/psb_zcgstab.f90 @@ -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 diff --git a/linsolve/impl/psb_zcgstabl.f90 b/linsolve/impl/psb_zcgstabl.f90 index aba579327..cf9ecbba3 100644 --- a/linsolve/impl/psb_zcgstabl.f90 +++ b/linsolve/impl/psb_zcgstabl.f90 @@ -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 diff --git a/linsolve/impl/psb_zfcg.F90 b/linsolve/impl/psb_zfcg.F90 index 41dbfe8df..63ac023be 100644 --- a/linsolve/impl/psb_zfcg.F90 +++ b/linsolve/impl/psb_zfcg.F90 @@ -163,9 +163,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_)& diff --git a/linsolve/impl/psb_zgcr.f90 b/linsolve/impl/psb_zgcr.f90 index 55feffd1f..9dc62b6f6 100644 --- a/linsolve/impl/psb_zgcr.f90 +++ b/linsolve/impl/psb_zgcr.f90 @@ -166,22 +166,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) diff --git a/linsolve/impl/psb_zkrylov.f90 b/linsolve/impl/psb_zkrylov.f90 index 091271d15..3baff8c34 100644 --- a/linsolve/impl/psb_zkrylov.f90 +++ b/linsolve/impl/psb_zkrylov.f90 @@ -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) diff --git a/linsolve/impl/psb_zminres.f90 b/linsolve/impl/psb_zminres.f90 new file mode 100644 index 000000000..4d9c27dbe --- /dev/null +++ b/linsolve/impl/psb_zminres.f90 @@ -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 diff --git a/linsolve/impl/psb_zrgmres.f90 b/linsolve/impl/psb_zrgmres.f90 index fc1450dcd..017704e30 100644 --- a/linsolve/impl/psb_zrgmres.f90 +++ b/linsolve/impl/psb_zrgmres.f90 @@ -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| @@ -129,7 +130,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, 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. @@ -171,24 +172,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 @@ -247,12 +259,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_ @@ -267,7 +282,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 @@ -318,27 +334,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) @@ -354,42 +375,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) ! + 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) rni = psb_geamax(w1,desc_a,info) @@ -397,8 +394,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 @@ -406,7 +402,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 @@ -414,28 +417,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) @@ -448,25 +441,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) @@ -496,6 +479,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 diff --git a/linsolve/psb_base_linsolve_conv_mod.f90 b/linsolve/psb_base_linsolve_conv_mod.f90 index 8e45614f4..0e45d11e9 100644 --- a/linsolve/psb_base_linsolve_conv_mod.f90 +++ b/linsolve/psb_base_linsolve_conv_mod.f90 @@ -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 diff --git a/test/pdegen/psb_d_pde2d.F90 b/test/pdegen/psb_d_pde2d.F90 index 183b8d160..907d0d3f8 100644 --- a/test/pdegen/psb_d_pde2d.F90 +++ b/test/pdegen/psb_d_pde2d.F90 @@ -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) diff --git a/test/pdegen/psb_d_pde3d.F90 b/test/pdegen/psb_d_pde3d.F90 index b0e692565..9a3ec5b68 100644 --- a/test/pdegen/psb_d_pde3d.F90 +++ b/test/pdegen/psb_d_pde3d.F90 @@ -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) diff --git a/test/pdegen/psb_s_pde2d.F90 b/test/pdegen/psb_s_pde2d.F90 index a4d724cab..fdef71e6b 100644 --- a/test/pdegen/psb_s_pde2d.F90 +++ b/test/pdegen/psb_s_pde2d.F90 @@ -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) diff --git a/test/pdegen/psb_s_pde3d.F90 b/test/pdegen/psb_s_pde3d.F90 index 87e631c72..ec4ee0f1f 100644 --- a/test/pdegen/psb_s_pde3d.F90 +++ b/test/pdegen/psb_s_pde3d.F90 @@ -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) diff --git a/util/psb_c_mmio_impl.f90 b/util/psb_c_mmio_impl.f90 index 68690d06c..ad653462a 100644 --- a/util/psb_c_mmio_impl.f90 +++ b/util/psb_c_mmio_impl.f90 @@ -36,12 +36,13 @@ subroutine mm_cvet_read(b, info, iunit, filename) integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename - integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j,infile - real(psb_spk_) :: bre, bim + integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, infile character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& & line*1024 + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -52,6 +53,7 @@ subroutine mm_cvet_read(b, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -76,15 +78,15 @@ subroutine mm_cvet_read(b, info, iunit, filename) read(line,fmt=*)nrow,ncol - if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'general')) then + if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'general')) then allocate(b(nrow),stat = ircode) if (ircode /= 0) goto 993 - do i=1, nrow - read(infile,fmt=*,end=902) bre,bim - b(i) = cmplx(bre,bim,kind=psb_spk_) + do i=1,nrow + read(infile,fmt=*,end=902) b(i) end do + end if ! read right hand sides - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -109,12 +111,13 @@ subroutine mm_cvet2_read(b, info, iunit, filename) integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename - integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j,infile - real(psb_spk_) :: bre, bim + integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, infile character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& & line*1024 + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -125,6 +128,7 @@ subroutine mm_cvet2_read(b, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -152,15 +156,10 @@ subroutine mm_cvet2_read(b, info, iunit, filename) if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'general')) then allocate(b(nrow,ncol),stat = ircode) if (ircode /= 0) goto 993 - do j=1, ncol - do i=1, nrow - read(infile,fmt=*,end=902) bre,bim - b(i,j) = cmplx(bre,bim,kind=psb_spk_) - end do - end do + read(infile,fmt=*,end=902) ((b(i,j), i=1,nrow),j=1,ncol) end if ! read right hand sides - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -189,8 +188,10 @@ subroutine mm_cvet2_write(b, header, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile character(len=80) :: frmtv + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then outfile=6 @@ -201,6 +202,7 @@ subroutine mm_cvet2_write(b, header, info, iunit, filename) outfile=99 endif open(outfile,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -209,17 +211,19 @@ subroutine mm_cvet2_write(b, header, info, iunit, filename) outfile=6 endif endif - - write(outfile,'(a)') '%%MatrixMarket matrix array complex general' + + write(outfile,'(a)') '%%MatrixMarket matrix array real general' write(outfile,'(a)') '% '//trim(header) write(outfile,'(a)') '% ' nrow = size(b,1) ncol = size(b,2) - write(outfile,*) nrow,ncol + write(outfile,*) nrow, ncol + + write(frmtv,'(a,i0,a)') '(',ncol,'(es26.18,1x))' - write(outfile,fmt='(2(es26.18,1x))') ((b(i,j), i=1,nrow),j=1,ncol) + write(outfile,fmt=frmtv) ((b(i,j), i=1,nrow),j=1,ncol) - if (outfile /= 6) close(outfile) + if (opened) close(outfile) return ! open failed @@ -241,8 +245,10 @@ subroutine mm_cvet1_write(b, header, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile character(len=80) :: frmtv + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then outfile=6 @@ -253,6 +259,7 @@ subroutine mm_cvet1_write(b, header, info, iunit, filename) outfile=99 endif open(outfile,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -262,20 +269,20 @@ subroutine mm_cvet1_write(b, header, info, iunit, filename) endif endif - write(outfile,'(a)') '%%MatrixMarket matrix array complex general' + write(outfile,'(a)') '%%MatrixMarket matrix array real general' write(outfile,'(a)') '% '//trim(header) write(outfile,'(a)') '% ' nrow = size(b,1) ncol = 1 write(outfile,*) nrow,ncol - write(frmtv,'(a,i0,a)') '(',2*ncol,'(es26.18,1x))' + write(frmtv,'(a,i0,a)') '(',ncol,'(es26.18,1x))' do i=1,size(b,1) write(outfile,frmtv) b(i) end do - if (outfile /= 6) close(outfile) + if (opened) close(outfile) return ! open failed @@ -298,7 +305,7 @@ subroutine mm_cvect_read(b, info, iunit, filename) complex(psb_spk_), allocatable :: bv(:) call mm_array_read(bv, info, iunit, filename) - call b%bld(bv) + if (info == 0) call b%bld(bv) end subroutine mm_cvect_read @@ -331,9 +338,10 @@ subroutine cmm_mat_read(a, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, nnzero integer(psb_ipk_) :: ircode, i,nzr,infile type(psb_c_coo_sparse_mat), allocatable :: acoo - real(psb_spk_) :: are, aim - info = psb_success_ + logical :: opened + info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -344,6 +352,7 @@ subroutine cmm_mat_read(a, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -369,24 +378,27 @@ subroutine cmm_mat_read(a, info, iunit, filename) allocate(acoo, stat=ircode) if (ircode /= 0) goto 993 - if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'general')) then + if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'general')) then call acoo%allocate(nrow,ncol,nnzero) do i=1,nnzero - read(infile,fmt=*,end=902) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_spk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i) end do call acoo%set_nzeros(nnzero) - call acoo%fix(info) - call a%mv_from(acoo) + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'general')) then + call acoo%allocate(nrow,ncol,nnzero) + do i=1,nnzero + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) + end do + acoo%val(:) = cone + call acoo%set_nzeros(nnzero) - else if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'symmetric')) then + else if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'symmetric')) then ! we are generally working with non-symmetric matrices, so ! we de-symmetrize what we are about to read call acoo%allocate(nrow,ncol,2*nnzero) do i=1,nnzero - read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_spk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i) end do nzr = nnzero do i=1,nnzero @@ -398,37 +410,34 @@ subroutine cmm_mat_read(a, info, iunit, filename) end if end do call acoo%set_nzeros(nzr) - call acoo%fix(info) - call a%mv_from(acoo) - - else if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'hermitian')) then - ! we are generally working with non-symmetric matrices, so - ! we de-symmetrize what we are about to read + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'symmetric')) then call acoo%allocate(nrow,ncol,2*nnzero) do i=1,nnzero - read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_spk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) end do + acoo%val(:) = cone nzr = nnzero do i=1,nnzero if (acoo%ia(i) /= acoo%ja(i)) then nzr = nzr + 1 - acoo%val(nzr) = conjg(acoo%val(i)) acoo%ia(nzr) = acoo%ja(i) acoo%ja(nzr) = acoo%ia(i) end if end do call acoo%set_nzeros(nzr) - call acoo%fix(info) - - call a%mv_from(acoo) else write(psb_err_unit,*) 'read_matrix: matrix type not yet supported' info=904 end if - if (infile /= 5) close(infile) + + if (info == 0) then + call acoo%fix(info) + call a%mv_from(acoo) + end if + + if (opened) close(infile) return ! open failed @@ -456,10 +465,11 @@ subroutine cmm_mat_write(a,mtitle,info,iunit,filename) integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename integer(psb_ipk_) :: iout + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then iout=6 @@ -470,6 +480,7 @@ subroutine cmm_mat_write(a,mtitle,info,iunit,filename) iout=99 endif open(iout,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -481,7 +492,7 @@ subroutine cmm_mat_write(a,mtitle,info,iunit,filename) call a%print(iout,head=mtitle) - if (iout /= 6) close(iout) + if (opened) close(iout) return @@ -492,6 +503,7 @@ subroutine cmm_mat_write(a,mtitle,info,iunit,filename) return end subroutine cmm_mat_write + subroutine lcmm_mat_read(a, info, iunit, filename) use psb_base_mod implicit none @@ -501,12 +513,13 @@ subroutine lcmm_mat_read(a, info, iunit, filename) character(len=*), optional, intent(in) :: filename character :: mmheader*15, fmt*15, object*10, type*10, sym*15 character(1024) :: line - integer(psb_lpk_) :: nrow, ncol, nnzero, i,nzr - integer(psb_ipk_) :: ircode,infile + integer(psb_lpk_) :: nrow, ncol, nnzero, i, nzr + integer(psb_ipk_) :: ircode, infile type(psb_lc_coo_sparse_mat), allocatable :: acoo - real(psb_spk_) :: are, aim - info = psb_success_ + logical :: opened + info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -517,6 +530,7 @@ subroutine lcmm_mat_read(a, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -542,24 +556,27 @@ subroutine lcmm_mat_read(a, info, iunit, filename) allocate(acoo, stat=ircode) if (ircode /= 0) goto 993 - if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'general')) then + if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'general')) then call acoo%allocate(nrow,ncol,nnzero) do i=1,nnzero - read(infile,fmt=*,end=902) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_spk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i) end do call acoo%set_nzeros(nnzero) - call acoo%fix(info) - call a%mv_from(acoo) + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'general')) then + call acoo%allocate(nrow,ncol,nnzero) + do i=1,nnzero + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) + end do + acoo%val(:) = cone + call acoo%set_nzeros(nnzero) - else if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'symmetric')) then + else if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'symmetric')) then ! we are generally working with non-symmetric matrices, so ! we de-symmetrize what we are about to read call acoo%allocate(nrow,ncol,2*nnzero) do i=1,nnzero - read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_spk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i) end do nzr = nnzero do i=1,nnzero @@ -571,37 +588,34 @@ subroutine lcmm_mat_read(a, info, iunit, filename) end if end do call acoo%set_nzeros(nzr) - call acoo%fix(info) - call a%mv_from(acoo) - - else if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'hermitian')) then - ! we are generally working with non-symmetric matrices, so - ! we de-symmetrize what we are about to read + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'symmetric')) then call acoo%allocate(nrow,ncol,2*nnzero) do i=1,nnzero - read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_spk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) end do + acoo%val(:) = cone nzr = nnzero do i=1,nnzero if (acoo%ia(i) /= acoo%ja(i)) then nzr = nzr + 1 - acoo%val(nzr) = conjg(acoo%val(i)) acoo%ia(nzr) = acoo%ja(i) acoo%ja(nzr) = acoo%ia(i) end if end do call acoo%set_nzeros(nzr) - call acoo%fix(info) - - call a%mv_from(acoo) else write(psb_err_unit,*) 'read_matrix: matrix type not yet supported' info=904 end if - if (infile /= 5) close(infile) + + if (info == 0) then + call acoo%fix(info) + call a%mv_from(acoo) + end if + + if (opened) close(infile) return ! open failed @@ -629,10 +643,10 @@ subroutine lcmm_mat_write(a,mtitle,info,iunit,filename) integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename integer(psb_ipk_) :: iout - + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then iout=6 @@ -643,6 +657,7 @@ subroutine lcmm_mat_write(a,mtitle,info,iunit,filename) iout=99 endif open(iout,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -654,7 +669,7 @@ subroutine lcmm_mat_write(a,mtitle,info,iunit,filename) call a%print(iout,head=mtitle) - if (iout /= 6) close(iout) + if (opened) close(iout) return diff --git a/util/psb_d_mmio_impl.f90 b/util/psb_d_mmio_impl.f90 index 374eeed92..4fc917a61 100644 --- a/util/psb_d_mmio_impl.f90 +++ b/util/psb_d_mmio_impl.f90 @@ -39,8 +39,10 @@ subroutine mm_dvet_read(b, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, infile character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& & line*1024 + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -51,6 +53,7 @@ subroutine mm_dvet_read(b, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -83,7 +86,7 @@ subroutine mm_dvet_read(b, info, iunit, filename) end do end if ! read right hand sides - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -111,8 +114,10 @@ subroutine mm_dvet2_read(b, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, infile character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& & line*1024 + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -123,6 +128,7 @@ subroutine mm_dvet2_read(b, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -153,7 +159,7 @@ subroutine mm_dvet2_read(b, info, iunit, filename) read(infile,fmt=*,end=902) ((b(i,j), i=1,nrow),j=1,ncol) end if ! read right hand sides - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -182,8 +188,10 @@ subroutine mm_dvet2_write(b, header, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile character(len=80) :: frmtv + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then outfile=6 @@ -194,6 +202,7 @@ subroutine mm_dvet2_write(b, header, info, iunit, filename) outfile=99 endif open(outfile,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -202,7 +211,7 @@ subroutine mm_dvet2_write(b, header, info, iunit, filename) outfile=6 endif endif - + write(outfile,'(a)') '%%MatrixMarket matrix array real general' write(outfile,'(a)') '% '//trim(header) write(outfile,'(a)') '% ' @@ -210,9 +219,11 @@ subroutine mm_dvet2_write(b, header, info, iunit, filename) ncol = size(b,2) write(outfile,*) nrow, ncol - write(outfile,fmt='(es26.18,1x)') ((b(i,j), i=1,nrow),j=1,ncol) + write(frmtv,'(a,i0,a)') '(',ncol,'(es26.18,1x))' + + write(outfile,fmt=frmtv) ((b(i,j), i=1,nrow),j=1,ncol) - if (outfile /= 6) close(outfile) + if (opened) close(outfile) return ! open failed @@ -234,8 +245,10 @@ subroutine mm_dvet1_write(b, header, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile character(len=80) :: frmtv + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then outfile=6 @@ -246,6 +259,7 @@ subroutine mm_dvet1_write(b, header, info, iunit, filename) outfile=99 endif open(outfile,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -268,7 +282,7 @@ subroutine mm_dvet1_write(b, header, info, iunit, filename) write(outfile,frmtv) b(i) end do - if (outfile /= 6) close(outfile) + if (opened) close(outfile) return ! open failed @@ -324,9 +338,10 @@ subroutine dmm_mat_read(a, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, nnzero integer(psb_ipk_) :: ircode, i,nzr,infile type(psb_d_coo_sparse_mat), allocatable :: acoo + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -337,6 +352,7 @@ subroutine dmm_mat_read(a, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -421,7 +437,7 @@ subroutine dmm_mat_read(a, info, iunit, filename) call a%mv_from(acoo) end if - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -449,10 +465,11 @@ subroutine dmm_mat_write(a,mtitle,info,iunit,filename) integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename integer(psb_ipk_) :: iout + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then iout=6 @@ -463,6 +480,7 @@ subroutine dmm_mat_write(a,mtitle,info,iunit,filename) iout=99 endif open(iout,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -474,7 +492,7 @@ subroutine dmm_mat_write(a,mtitle,info,iunit,filename) call a%print(iout,head=mtitle) - if (iout /= 6) close(iout) + if (opened) close(iout) return @@ -498,9 +516,10 @@ subroutine ldmm_mat_read(a, info, iunit, filename) integer(psb_lpk_) :: nrow, ncol, nnzero, i, nzr integer(psb_ipk_) :: ircode, infile type(psb_ld_coo_sparse_mat), allocatable :: acoo + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -511,6 +530,7 @@ subroutine ldmm_mat_read(a, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -595,7 +615,7 @@ subroutine ldmm_mat_read(a, info, iunit, filename) call a%mv_from(acoo) end if - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -623,10 +643,10 @@ subroutine ldmm_mat_write(a,mtitle,info,iunit,filename) integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename integer(psb_ipk_) :: iout - + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then iout=6 @@ -637,6 +657,7 @@ subroutine ldmm_mat_write(a,mtitle,info,iunit,filename) iout=99 endif open(iout,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -648,7 +669,7 @@ subroutine ldmm_mat_write(a,mtitle,info,iunit,filename) call a%print(iout,head=mtitle) - if (iout /= 6) close(iout) + if (opened) close(iout) return @@ -658,5 +679,3 @@ subroutine ldmm_mat_write(a,mtitle,info,iunit,filename) write(psb_err_unit,*) 'Error while opening ',filename return end subroutine ldmm_mat_write - - diff --git a/util/psb_i_mmio_impl.f90 b/util/psb_i_mmio_impl.f90 new file mode 100644 index 000000000..64d02bcd7 --- /dev/null +++ b/util/psb_i_mmio_impl.f90 @@ -0,0 +1,328 @@ +! +! 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. +! +! +subroutine mm_ivet_read(b, info, iunit, filename) + use psb_base_mod + implicit none + integer(psb_ipk_), allocatable, intent(out) :: b(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: iunit + character(len=*), optional, intent(in) :: filename + integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, infile + character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& + & line*1024 + logical :: opened + + info = psb_success_ + opened = .false. + if (present(filename)) then + if (filename == '-') then + infile=5 + else + if (present(iunit)) then + infile=iunit + else + infile=99 + endif + open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. + endif + else + if (present(iunit)) then + infile=iunit + else + infile=5 + endif + endif + + read(infile,fmt=*, end=902) mmheader, object, fmt, type, sym + + if ( (object /= 'matrix').or.(fmt /= 'array')) then + write(psb_err_unit,*) 'read_rhs: input file type not yet supported' + info = -3 + return + end if + + do + read(infile,fmt='(a)') line + if (line(1:1) /= '%') exit + end do + + read(line,fmt=*)nrow,ncol + + if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'general')) then + allocate(b(nrow),stat = ircode) + if (ircode /= 0) goto 993 + do i=1,nrow + read(infile,fmt=*,end=902) b(i) + end do + + end if ! read right hand sides + if (opened) close(infile) + + return + ! open failed +901 write(psb_err_unit,*) 'mm_vet_read: could not open file ',& + & infile,' for input' + info = -1 + return + +902 write(psb_err_unit,*) 'mmv_vet_read: unexpected end of file ',infile,& + & ' during input' + info = -2 + return +993 write(psb_err_unit,*) 'mm_vet_read: memory allocation failure' + info = -3 + return +end subroutine mm_ivet_read + +subroutine mm_ivet2_read(b, info, iunit, filename) + use psb_base_mod + implicit none + integer(psb_ipk_), allocatable, intent(out) :: b(:,:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: iunit + character(len=*), optional, intent(in) :: filename + integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, infile + character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& + & line*1024 + logical :: opened + + info = psb_success_ + opened = .false. + if (present(filename)) then + if (filename == '-') then + infile=5 + else + if (present(iunit)) then + infile=iunit + else + infile=99 + endif + open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. + endif + else + if (present(iunit)) then + infile=iunit + else + infile=5 + endif + endif + + read(infile,fmt=*, end=902) mmheader, object, fmt, type, sym + + if ( (object /= 'matrix').or.(fmt /= 'array')) then + write(psb_err_unit,*) 'read_rhs: input file type not yet supported' + info = -3 + return + end if + + do + read(infile,fmt='(a)') line + if (line(1:1) /= '%') exit + end do + + read(line,fmt=*)nrow,ncol + + if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'general')) then + allocate(b(nrow,ncol),stat = ircode) + if (ircode /= 0) goto 993 + read(infile,fmt=*,end=902) ((b(i,j), i=1,nrow),j=1,ncol) + + end if ! read right hand sides + if (opened) close(infile) + + return + ! open failed +901 write(psb_err_unit,*) 'mm_vet_read: could not open file ',& + & infile,' for input' + info = -1 + return + +902 write(psb_err_unit,*) 'mmv_vet_read: unexpected end of file ',infile,& + & ' during input' + info = -2 + return +993 write(psb_err_unit,*) 'mm_vet_read: memory allocation failure' + info = -3 + return +end subroutine mm_ivet2_read + +subroutine mm_ivet2_write(b, header, info, iunit, filename) + use psb_base_mod + implicit none + integer(psb_ipk_), intent(in) :: b(:,:) + character(len=*), intent(in) :: header + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: iunit + character(len=*), optional, intent(in) :: filename + integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile + + character(len=80) :: frmtv + logical :: opened + + info = psb_success_ + opened = .false. + if (present(filename)) then + if (filename == '-') then + outfile=6 + else + if (present(iunit)) then + outfile=iunit + else + outfile=99 + endif + open(outfile,file=filename, err=901, action='WRITE') + opened = .true. + endif + else + if (present(iunit)) then + outfile=iunit + else + outfile=6 + endif + endif + + write(outfile,'(a)') '%%MatrixMarket matrix array real general' + write(outfile,'(a)') '% '//trim(header) + write(outfile,'(a)') '% ' + nrow = size(b,1) + ncol = size(b,2) + write(outfile,*) nrow, ncol + + write(frmtv,'(a,i0,a)') '(',ncol,'(es26.18,1x))' + + write(outfile,fmt=frmtv) ((b(i,j), i=1,nrow),j=1,ncol) + + if (opened) close(outfile) + + return + ! open failed +901 write(psb_err_unit,*) 'mm_vet_write: could not open file ',& + & outfile,' for output' + info = -1 + return + +end subroutine mm_ivet2_write + +subroutine mm_ivet1_write(b, header, info, iunit, filename) + use psb_base_mod + implicit none + integer(psb_ipk_), intent(in) :: b(:) + character(len=*), intent(in) :: header + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: iunit + character(len=*), optional, intent(in) :: filename + integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile + + character(len=80) :: frmtv + logical :: opened + + info = psb_success_ + opened = .false. + if (present(filename)) then + if (filename == '-') then + outfile=6 + else + if (present(iunit)) then + outfile=iunit + else + outfile=99 + endif + open(outfile,file=filename, err=901, action='WRITE') + opened = .true. + endif + else + if (present(iunit)) then + outfile=iunit + else + outfile=6 + endif + endif + + write(outfile,'(a)') '%%MatrixMarket matrix array real general' + write(outfile,'(a)') '% '//trim(header) + write(outfile,'(a)') '% ' + nrow = size(b,1) + ncol = 1 + write(outfile,*) nrow,ncol + + write(frmtv,'(a,i0,a)') '(',ncol,'(es26.18,1x))' + + do i=1,size(b,1) + write(outfile,frmtv) b(i) + end do + + if (opened) close(outfile) + + return + ! open failed +901 write(psb_err_unit,*) 'mm_vet_write: could not open file ',& + & outfile,' for output' + info = -1 + return + +end subroutine mm_ivet1_write + +subroutine mm_ivect_read(b, info, iunit, filename) + use psb_base_mod + use psb_mmio_mod, psb_protect_name => mm_ivect_read + implicit none + type(psb_i_vect_type), intent(inout) :: b + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: iunit + character(len=*), optional, intent(in) :: filename + ! + integer(psb_ipk_), allocatable :: bv(:) + + call mm_array_read(bv, info, iunit, filename) + if (info == 0) call b%bld(bv) + +end subroutine mm_ivect_read + +subroutine mm_ivect_write(b, header, info, iunit, filename) + use psb_base_mod + use psb_mmio_mod, psb_protect_name => mm_ivect_write + implicit none + type(psb_i_vect_type), intent(inout) :: b + character(len=*), intent(in) :: header + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: iunit + character(len=*), optional, intent(in) :: filename + info = psb_success_ + if (.not.allocated(b%v)) return + call b%sync() + + call mm_array_write(b%v%v,header,info,iunit,filename) + +end subroutine mm_ivect_write + diff --git a/util/psb_s_mmio_impl.f90 b/util/psb_s_mmio_impl.f90 index 94e886849..0ffd9a93d 100644 --- a/util/psb_s_mmio_impl.f90 +++ b/util/psb_s_mmio_impl.f90 @@ -36,11 +36,13 @@ subroutine mm_svet_read(b, info, iunit, filename) integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename - integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j,infile + integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, infile character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& & line*1024 + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -51,6 +53,7 @@ subroutine mm_svet_read(b, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -83,8 +86,7 @@ subroutine mm_svet_read(b, info, iunit, filename) end do end if ! read right hand sides - - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -109,11 +111,13 @@ subroutine mm_svet2_read(b, info, iunit, filename) integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename - integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j,infile + integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, infile character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& & line*1024 + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -124,6 +128,7 @@ subroutine mm_svet2_read(b, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -154,8 +159,7 @@ subroutine mm_svet2_read(b, info, iunit, filename) read(infile,fmt=*,end=902) ((b(i,j), i=1,nrow),j=1,ncol) end if ! read right hand sides - - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -184,8 +188,10 @@ subroutine mm_svet2_write(b, header, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile character(len=80) :: frmtv + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then outfile=6 @@ -196,6 +202,7 @@ subroutine mm_svet2_write(b, header, info, iunit, filename) outfile=99 endif open(outfile,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -204,17 +211,19 @@ subroutine mm_svet2_write(b, header, info, iunit, filename) outfile=6 endif endif - + write(outfile,'(a)') '%%MatrixMarket matrix array real general' write(outfile,'(a)') '% '//trim(header) write(outfile,'(a)') '% ' nrow = size(b,1) ncol = size(b,2) - write(outfile,*) nrow,ncol + write(outfile,*) nrow, ncol + + write(frmtv,'(a,i0,a)') '(',ncol,'(es26.18,1x))' - write(outfile,fmt='(es26.18,1x)') ((b(i,j), i=1,nrow),j=1,ncol) + write(outfile,fmt=frmtv) ((b(i,j), i=1,nrow),j=1,ncol) - if (outfile /= 6) close(outfile) + if (opened) close(outfile) return ! open failed @@ -236,8 +245,10 @@ subroutine mm_svet1_write(b, header, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile character(len=80) :: frmtv + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then outfile=6 @@ -248,6 +259,7 @@ subroutine mm_svet1_write(b, header, info, iunit, filename) outfile=99 endif open(outfile,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -270,7 +282,7 @@ subroutine mm_svet1_write(b, header, info, iunit, filename) write(outfile,frmtv) b(i) end do - if (outfile /= 6) close(outfile) + if (opened) close(outfile) return ! open failed @@ -326,9 +338,10 @@ subroutine smm_mat_read(a, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, nnzero integer(psb_ipk_) :: ircode, i,nzr,infile type(psb_s_coo_sparse_mat), allocatable :: acoo + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -339,6 +352,7 @@ subroutine smm_mat_read(a, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -370,9 +384,14 @@ subroutine smm_mat_read(a, info, iunit, filename) read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i) end do call acoo%set_nzeros(nnzero) - call acoo%fix(info) - call a%mv_from(acoo) + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'general')) then + call acoo%allocate(nrow,ncol,nnzero) + do i=1,nnzero + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) + end do + acoo%val(:) = sone + call acoo%set_nzeros(nnzero) else if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'symmetric')) then ! we are generally working with non-symmetric matrices, so @@ -391,17 +410,34 @@ subroutine smm_mat_read(a, info, iunit, filename) end if end do call acoo%set_nzeros(nzr) - call acoo%fix(info) - call a%mv_from(acoo) + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'symmetric')) then + call acoo%allocate(nrow,ncol,2*nnzero) + do i=1,nnzero + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) + end do + acoo%val(:) = sone + nzr = nnzero + do i=1,nnzero + if (acoo%ia(i) /= acoo%ja(i)) then + nzr = nzr + 1 + acoo%ia(nzr) = acoo%ja(i) + acoo%ja(nzr) = acoo%ia(i) + end if + end do + call acoo%set_nzeros(nzr) else write(psb_err_unit,*) 'read_matrix: matrix type not yet supported' info=904 end if + if (info == 0) then + call acoo%fix(info) + call a%mv_from(acoo) + end if - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -429,10 +465,11 @@ subroutine smm_mat_write(a,mtitle,info,iunit,filename) integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename integer(psb_ipk_) :: iout + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then iout=6 @@ -443,6 +480,7 @@ subroutine smm_mat_write(a,mtitle,info,iunit,filename) iout=99 endif open(iout,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -454,7 +492,7 @@ subroutine smm_mat_write(a,mtitle,info,iunit,filename) call a%print(iout,head=mtitle) - if (iout /= 6) close(iout) + if (opened) close(iout) return @@ -465,6 +503,7 @@ subroutine smm_mat_write(a,mtitle,info,iunit,filename) return end subroutine smm_mat_write + subroutine lsmm_mat_read(a, info, iunit, filename) use psb_base_mod implicit none @@ -474,12 +513,13 @@ subroutine lsmm_mat_read(a, info, iunit, filename) character(len=*), optional, intent(in) :: filename character :: mmheader*15, fmt*15, object*10, type*10, sym*15 character(1024) :: line - integer(psb_lpk_) :: nrow, ncol, nnzero, i,nzr - integer(psb_ipk_) :: ircode,infile + integer(psb_lpk_) :: nrow, ncol, nnzero, i, nzr + integer(psb_ipk_) :: ircode, infile type(psb_ls_coo_sparse_mat), allocatable :: acoo + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -490,6 +530,7 @@ subroutine lsmm_mat_read(a, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -521,9 +562,14 @@ subroutine lsmm_mat_read(a, info, iunit, filename) read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i) end do call acoo%set_nzeros(nnzero) - call acoo%fix(info) - call a%mv_from(acoo) + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'general')) then + call acoo%allocate(nrow,ncol,nnzero) + do i=1,nnzero + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) + end do + acoo%val(:) = sone + call acoo%set_nzeros(nnzero) else if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'symmetric')) then ! we are generally working with non-symmetric matrices, so @@ -542,17 +588,34 @@ subroutine lsmm_mat_read(a, info, iunit, filename) end if end do call acoo%set_nzeros(nzr) - call acoo%fix(info) - call a%mv_from(acoo) + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'symmetric')) then + call acoo%allocate(nrow,ncol,2*nnzero) + do i=1,nnzero + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) + end do + acoo%val(:) = sone + nzr = nnzero + do i=1,nnzero + if (acoo%ia(i) /= acoo%ja(i)) then + nzr = nzr + 1 + acoo%ia(nzr) = acoo%ja(i) + acoo%ja(nzr) = acoo%ia(i) + end if + end do + call acoo%set_nzeros(nzr) else write(psb_err_unit,*) 'read_matrix: matrix type not yet supported' info=904 end if + if (info == 0) then + call acoo%fix(info) + call a%mv_from(acoo) + end if - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -580,10 +643,10 @@ subroutine lsmm_mat_write(a,mtitle,info,iunit,filename) integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename integer(psb_ipk_) :: iout - + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then iout=6 @@ -594,6 +657,7 @@ subroutine lsmm_mat_write(a,mtitle,info,iunit,filename) iout=99 endif open(iout,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -605,7 +669,7 @@ subroutine lsmm_mat_write(a,mtitle,info,iunit,filename) call a%print(iout,head=mtitle) - if (iout /= 6) close(iout) + if (opened) close(iout) return diff --git a/util/psb_z_mmio_impl.f90 b/util/psb_z_mmio_impl.f90 index 948de2833..aa66f0edd 100644 --- a/util/psb_z_mmio_impl.f90 +++ b/util/psb_z_mmio_impl.f90 @@ -36,12 +36,13 @@ subroutine mm_zvet_read(b, info, iunit, filename) integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename - integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j,infile - real(psb_dpk_) :: bre, bim + integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, infile character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& & line*1024 + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -52,6 +53,7 @@ subroutine mm_zvet_read(b, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -76,15 +78,15 @@ subroutine mm_zvet_read(b, info, iunit, filename) read(line,fmt=*)nrow,ncol - if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'general')) then + if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'general')) then allocate(b(nrow),stat = ircode) if (ircode /= 0) goto 993 - do i=1, nrow - read(infile,fmt=*,end=902) bre,bim - b(i) = cmplx(bre,bim,kind=psb_dpk_) + do i=1,nrow + read(infile,fmt=*,end=902) b(i) end do + end if ! read right hand sides - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -109,12 +111,13 @@ subroutine mm_zvet2_read(b, info, iunit, filename) integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename - integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j,infile - real(psb_dpk_) :: bre, bim + integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, infile character :: mmheader*15, fmt*15, object*10, type*10, sym*15,& & line*1024 + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -125,6 +128,7 @@ subroutine mm_zvet2_read(b, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -149,18 +153,13 @@ subroutine mm_zvet2_read(b, info, iunit, filename) read(line,fmt=*)nrow,ncol - if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'general')) then + if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'general')) then allocate(b(nrow,ncol),stat = ircode) if (ircode /= 0) goto 993 - do j=1, ncol - do i=1, nrow - read(infile,fmt=*,end=902) bre,bim - b(i,j) = cmplx(bre,bim,kind=psb_dpk_) - end do - end do + read(infile,fmt=*,end=902) ((b(i,j), i=1,nrow),j=1,ncol) end if ! read right hand sides - if (infile /= 5) close(infile) + if (opened) close(infile) return ! open failed @@ -189,8 +188,10 @@ subroutine mm_zvet2_write(b, header, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile character(len=80) :: frmtv + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then outfile=6 @@ -201,6 +202,7 @@ subroutine mm_zvet2_write(b, header, info, iunit, filename) outfile=99 endif open(outfile,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -209,17 +211,19 @@ subroutine mm_zvet2_write(b, header, info, iunit, filename) outfile=6 endif endif - - write(outfile,'(a)') '%%MatrixMarket matrix array complex general' + + write(outfile,'(a)') '%%MatrixMarket matrix array real general' write(outfile,'(a)') '% '//trim(header) write(outfile,'(a)') '% ' nrow = size(b,1) ncol = size(b,2) - write(outfile,*) nrow,ncol + write(outfile,*) nrow, ncol + + write(frmtv,'(a,i0,a)') '(',ncol,'(es26.18,1x))' - write(outfile,fmt='(2(es26.18,1x))') ((b(i,j), i=1,nrow),j=1,ncol) + write(outfile,fmt=frmtv) ((b(i,j), i=1,nrow),j=1,ncol) - if (outfile /= 6) close(outfile) + if (opened) close(outfile) return ! open failed @@ -241,8 +245,10 @@ subroutine mm_zvet1_write(b, header, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, i,root, np, me, ircode, j, outfile character(len=80) :: frmtv + logical :: opened info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then outfile=6 @@ -253,6 +259,7 @@ subroutine mm_zvet1_write(b, header, info, iunit, filename) outfile=99 endif open(outfile,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -269,13 +276,13 @@ subroutine mm_zvet1_write(b, header, info, iunit, filename) ncol = 1 write(outfile,*) nrow,ncol - write(frmtv,'(a,i0,a)') '(',2*ncol,'(es26.18,1x))' + write(frmtv,'(a,i0,a)') '(',ncol,'(es26.18,1x))' do i=1,size(b,1) write(outfile,frmtv) b(i) end do - if (outfile /= 6) close(outfile) + if (opened) close(outfile) return ! open failed @@ -331,9 +338,10 @@ subroutine zmm_mat_read(a, info, iunit, filename) integer(psb_ipk_) :: nrow, ncol, nnzero integer(psb_ipk_) :: ircode, i,nzr,infile type(psb_z_coo_sparse_mat), allocatable :: acoo - real(psb_dpk_) :: are, aim - info = psb_success_ + logical :: opened + info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -344,6 +352,7 @@ subroutine zmm_mat_read(a, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -369,24 +378,27 @@ subroutine zmm_mat_read(a, info, iunit, filename) allocate(acoo, stat=ircode) if (ircode /= 0) goto 993 - if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'general')) then + if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'general')) then call acoo%allocate(nrow,ncol,nnzero) do i=1,nnzero - read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_dpk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i) end do call acoo%set_nzeros(nnzero) - call acoo%fix(info) - call a%mv_from(acoo) + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'general')) then + call acoo%allocate(nrow,ncol,nnzero) + do i=1,nnzero + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) + end do + acoo%val(:) = zone + call acoo%set_nzeros(nnzero) - else if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'symmetric')) then + else if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'symmetric')) then ! we are generally working with non-symmetric matrices, so ! we de-symmetrize what we are about to read call acoo%allocate(nrow,ncol,2*nnzero) do i=1,nnzero - read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_dpk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i) end do nzr = nnzero do i=1,nnzero @@ -398,37 +410,34 @@ subroutine zmm_mat_read(a, info, iunit, filename) end if end do call acoo%set_nzeros(nzr) - call acoo%fix(info) - - call a%mv_from(acoo) - else if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'hermitian')) then - ! we are generally working with non-symmetric matrices, so - ! we de-symmetrize what we are about to read + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'symmetric')) then call acoo%allocate(nrow,ncol,2*nnzero) do i=1,nnzero - read(infile,fmt=*,end=902) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_dpk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) end do + acoo%val(:) = zone nzr = nnzero do i=1,nnzero if (acoo%ia(i) /= acoo%ja(i)) then nzr = nzr + 1 - acoo%val(nzr) = conjg(acoo%val(i)) acoo%ia(nzr) = acoo%ja(i) acoo%ja(nzr) = acoo%ia(i) end if end do call acoo%set_nzeros(nzr) - call acoo%fix(info) - - call a%mv_from(acoo) else write(psb_err_unit,*) 'read_matrix: matrix type not yet supported' info=904 end if - if (infile /= 5) close(infile) + + if (info == 0) then + call acoo%fix(info) + call a%mv_from(acoo) + end if + + if (opened) close(infile) return ! open failed @@ -446,6 +455,7 @@ subroutine zmm_mat_read(a, info, iunit, filename) return end subroutine zmm_mat_read + subroutine zmm_mat_write(a,mtitle,info,iunit,filename) use psb_base_mod implicit none @@ -455,10 +465,11 @@ subroutine zmm_mat_write(a,mtitle,info,iunit,filename) integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename integer(psb_ipk_) :: iout + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then iout=6 @@ -469,6 +480,7 @@ subroutine zmm_mat_write(a,mtitle,info,iunit,filename) iout=99 endif open(iout,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -480,7 +492,7 @@ subroutine zmm_mat_write(a,mtitle,info,iunit,filename) call a%print(iout,head=mtitle) - if (iout /= 6) close(iout) + if (opened) close(iout) return @@ -501,12 +513,13 @@ subroutine lzmm_mat_read(a, info, iunit, filename) character(len=*), optional, intent(in) :: filename character :: mmheader*15, fmt*15, object*10, type*10, sym*15 character(1024) :: line - integer(psb_lpk_) :: nrow, ncol, nnzero, i,nzr - integer(psb_ipk_) :: ircode,infile + integer(psb_lpk_) :: nrow, ncol, nnzero, i, nzr + integer(psb_ipk_) :: ircode, infile type(psb_lz_coo_sparse_mat), allocatable :: acoo - real(psb_dpk_) :: are, aim - info = psb_success_ + logical :: opened + info = psb_success_ + opened = .false. if (present(filename)) then if (filename == '-') then infile=5 @@ -517,6 +530,7 @@ subroutine lzmm_mat_read(a, info, iunit, filename) infile=99 endif open(infile,file=filename, status='OLD', err=901, action='READ') + opened = .true. endif else if (present(iunit)) then @@ -542,24 +556,27 @@ subroutine lzmm_mat_read(a, info, iunit, filename) allocate(acoo, stat=ircode) if (ircode /= 0) goto 993 - if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'general')) then + if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'general')) then call acoo%allocate(nrow,ncol,nnzero) do i=1,nnzero - read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_dpk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i) end do call acoo%set_nzeros(nnzero) - call acoo%fix(info) - call a%mv_from(acoo) + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'general')) then + call acoo%allocate(nrow,ncol,nnzero) + do i=1,nnzero + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) + end do + acoo%val(:) = zone + call acoo%set_nzeros(nnzero) - else if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'symmetric')) then + else if ((psb_tolower(type) == 'real').and.(psb_tolower(sym) == 'symmetric')) then ! we are generally working with non-symmetric matrices, so ! we de-symmetrize what we are about to read call acoo%allocate(nrow,ncol,2*nnzero) do i=1,nnzero - read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_dpk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i),acoo%val(i) end do nzr = nnzero do i=1,nnzero @@ -571,37 +588,34 @@ subroutine lzmm_mat_read(a, info, iunit, filename) end if end do call acoo%set_nzeros(nzr) - call acoo%fix(info) - call a%mv_from(acoo) - - else if ((psb_tolower(type) == 'complex').and.(psb_tolower(sym) == 'hermitian')) then - ! we are generally working with non-symmetric matrices, so - ! we de-symmetrize what we are about to read + else if ((psb_tolower(type) == 'pattern').and.(psb_tolower(sym) == 'symmetric')) then call acoo%allocate(nrow,ncol,2*nnzero) do i=1,nnzero - read(infile,fmt=*,end=902) acoo%ia(i),acoo%ja(i),are,aim - acoo%val(i) = cmplx(are,aim,kind=psb_dpk_) + read(infile,fmt=*,end=902,err=905) acoo%ia(i),acoo%ja(i) end do + acoo%val(:) = zone nzr = nnzero do i=1,nnzero if (acoo%ia(i) /= acoo%ja(i)) then nzr = nzr + 1 - acoo%val(nzr) = conjg(acoo%val(i)) acoo%ia(nzr) = acoo%ja(i) acoo%ja(nzr) = acoo%ia(i) end if end do call acoo%set_nzeros(nzr) - call acoo%fix(info) - - call a%mv_from(acoo) else write(psb_err_unit,*) 'read_matrix: matrix type not yet supported' info=904 end if - if (infile /= 5) close(infile) + + if (info == 0) then + call acoo%fix(info) + call a%mv_from(acoo) + end if + + if (opened) close(infile) return ! open failed @@ -619,6 +633,7 @@ subroutine lzmm_mat_read(a, info, iunit, filename) return end subroutine lzmm_mat_read + subroutine lzmm_mat_write(a,mtitle,info,iunit,filename) use psb_base_mod implicit none @@ -628,10 +643,10 @@ subroutine lzmm_mat_write(a,mtitle,info,iunit,filename) integer(psb_ipk_), optional, intent(in) :: iunit character(len=*), optional, intent(in) :: filename integer(psb_ipk_) :: iout - + logical :: opened info = psb_success_ - + opened = .false. if (present(filename)) then if (filename == '-') then iout=6 @@ -642,6 +657,7 @@ subroutine lzmm_mat_write(a,mtitle,info,iunit,filename) iout=99 endif open(iout,file=filename, err=901, action='WRITE') + opened = .true. endif else if (present(iunit)) then @@ -653,7 +669,7 @@ subroutine lzmm_mat_write(a,mtitle,info,iunit,filename) call a%print(iout,head=mtitle) - if (iout /= 6) close(iout) + if (opened) close(iout) return @@ -663,4 +679,3 @@ subroutine lzmm_mat_write(a,mtitle,info,iunit,filename) write(psb_err_unit,*) 'Error while opening ',filename return end subroutine lzmm_mat_write -