diff --git a/mlprec/impl/Makefile b/mlprec/impl/Makefile index 87465261..b17b5449 100644 --- a/mlprec/impl/Makefile +++ b/mlprec/impl/Makefile @@ -25,26 +25,29 @@ DINNEROBJS= mld_dcoarse_bld.o mld_dmlprec_bld.o \ mld_dilu0_fact.o mld_diluk_fact.o mld_dilut_fact.o mld_daggrmap_bld.o \ mld_d_dec_map_bld.o mld_dmlprec_aply.o mld_daggrmat_asb.o \ $(DMPFOBJS) mld_d_base_solver_impl.o mld_d_base_smoother_impl.o mld_d_onelev_impl.o\ - mld_d_as_smoother_impl.o mld_d_jac_smoother_impl.o + mld_d_as_smoother_impl.o mld_d_jac_smoother_impl.o \ + mld_d_diag_solver_impl.o mld_d_id_solver_impl.o SINNEROBJS= mld_scoarse_bld.o mld_smlprec_bld.o \ mld_silu0_fact.o mld_siluk_fact.o mld_silut_fact.o mld_saggrmap_bld.o \ mld_s_dec_map_bld.o mld_smlprec_aply.o mld_saggrmat_asb.o \ $(SMPFOBJS) mld_s_base_solver_impl.o mld_s_base_smoother_impl.o mld_s_onelev_impl.o\ - mld_s_as_smoother_impl.o mld_s_jac_smoother_impl.o + mld_s_as_smoother_impl.o mld_s_jac_smoother_impl.o \ + mld_s_diag_solver_impl.o mld_s_id_solver_impl.o ZINNEROBJS= mld_zcoarse_bld.o mld_zmlprec_bld.o \ mld_zilu0_fact.o mld_ziluk_fact.o mld_zilut_fact.o mld_zaggrmap_bld.o \ mld_z_dec_map_bld.o mld_zmlprec_aply.o mld_zaggrmat_asb.o \ $(ZMPFOBJS) mld_z_base_solver_impl.o mld_z_base_smoother_impl.o mld_z_onelev_impl.o\ - mld_z_as_smoother_impl.o mld_z_jac_smoother_impl.o + mld_z_as_smoother_impl.o mld_z_jac_smoother_impl.o \ + mld_z_diag_solver_impl.o mld_z_id_solver_impl.o CINNEROBJS= mld_ccoarse_bld.o mld_cmlprec_bld.o \ mld_cilu0_fact.o mld_ciluk_fact.o mld_cilut_fact.o mld_caggrmap_bld.o \ mld_c_dec_map_bld.o mld_cmlprec_aply.o mld_caggrmat_asb.o \ $(CMPFOBJS) mld_c_base_solver_impl.o mld_c_base_smoother_impl.o mld_c_onelev_impl.o\ - mld_c_as_smoother_impl.o mld_c_jac_smoother_impl.o - + mld_c_as_smoother_impl.o mld_c_jac_smoother_impl.o \ + mld_c_diag_solver_impl.o mld_c_id_solver_impl.o INNEROBJS= $(SINNEROBJS) $(DINNEROBJS) $(CINNEROBJS) $(ZINNEROBJS) diff --git a/mlprec/impl/mld_c_diag_solver_impl.f90 b/mlprec/impl/mld_c_diag_solver_impl.f90 new file mode 100644 index 00000000..0cbbc2ad --- /dev/null +++ b/mlprec/impl/mld_c_diag_solver_impl.f90 @@ -0,0 +1,413 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.0) +!!$ +!!$ (C) Copyright 2008,2009,2010, 2010 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ 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 MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific 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 MLD2P4 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 mld_c_diag_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) + use psb_base_mod + use mld_c_diag_solver, mld_protect_name => mld_c_diag_solver_apply_vect + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_c_diag_solver_type), intent(inout) :: sv + type(psb_c_vect_type), intent(inout) :: x + type(psb_c_vect_type), intent(inout) :: y + complex(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_spk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + + integer :: n_row,n_col + complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) + integer :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='c_diag_solver_apply' + + call psb_erractionsave(err_act) + + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T','C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + n_row = desc_data%get_local_rows() + n_col = desc_data%get_local_cols() + if (x%get_nrows() < n_row) then + info = 36 + call psb_errpush(info,name,i_err=(/2,n_row,0,0,0/)) + goto 9999 + end if + if (y%get_nrows() < n_row) then + info = 36 + call psb_errpush(info,name,i_err=(/3,n_row,0,0,0/)) + goto 9999 + end if + if (.not.allocated(sv%dv)) then + info = 1124 + call psb_errpush(info,name,a_err="preconditioner: D") + goto 9999 + end if + if (sv%dv%get_nrows() < n_row) then + info = 1124 + call psb_errpush(info,name,a_err="preconditioner: D") + goto 9999 + end if + + call y%mlt(alpha,sv%dv,x,beta,info,conjgx=trans_) + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='vect%mlt') + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + +end subroutine mld_c_diag_solver_apply_vect + +subroutine mld_c_diag_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) + use psb_base_mod + use mld_c_diag_solver, mld_protect_name => mld_c_diag_solver_apply + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_c_diag_solver_type), intent(in) :: sv + complex(psb_spk_), intent(inout) :: x(:) + complex(psb_spk_), intent(inout) :: y(:) + complex(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_spk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + + integer :: n_row,n_col + complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) + integer :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='c_diag_solver_apply' + + call psb_erractionsave(err_act) + + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T','C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + n_row = desc_data%get_local_rows() + n_col = desc_data%get_local_cols() + + if (trans_ == 'C') then + if (beta == czero) then + + if (alpha == czero) then + y(1:n_row) = czero + else if (alpha == cone) then + do i=1, n_row + y(i) = conjg(sv%d(i)) * x(i) + end do + else if (alpha == -cone) then + do i=1, n_row + y(i) = -conjg(sv%d(i)) * x(i) + end do + else + do i=1, n_row + y(i) = alpha * conjg(sv%d(i)) * x(i) + end do + end if + + else if (beta == cone) then + + if (alpha == czero) then + !y(1:n_row) = czero + else if (alpha == cone) then + do i=1, n_row + y(i) = conjg(sv%d(i)) * x(i) + y(i) + end do + else if (alpha == -cone) then + do i=1, n_row + y(i) = -conjg(sv%d(i)) * x(i) + y(i) + end do + else + do i=1, n_row + y(i) = alpha * conjg(sv%d(i)) * x(i) + y(i) + end do + end if + + else if (beta == -cone) then + + if (alpha == czero) then + y(1:n_row) = -y(1:n_row) + else if (alpha == cone) then + do i=1, n_row + y(i) = conjg(sv%d(i)) * x(i) - y(i) + end do + else if (alpha == -cone) then + do i=1, n_row + y(i) = -conjg(sv%d(i)) * x(i) - y(i) + end do + else + do i=1, n_row + y(i) = alpha * conjg(sv%d(i)) * x(i) - y(i) + end do + end if + + else + + if (alpha == czero) then + y(1:n_row) = beta *y(1:n_row) + else if (alpha == cone) then + do i=1, n_row + y(i) = conjg(sv%d(i)) * x(i) + beta*y(i) + end do + else if (alpha == -cone) then + do i=1, n_row + y(i) = -conjg(sv%d(i)) * x(i) + beta*y(i) + end do + else + do i=1, n_row + y(i) = alpha * conjg(sv%d(i)) * x(i) + beta*y(i) + end do + end if + + end if + + else if (trans_ /= 'C') then + + if (beta == czero) then + + if (alpha == czero) then + y(1:n_row) = czero + else if (alpha == cone) then + do i=1, n_row + y(i) = sv%d(i) * x(i) + end do + else if (alpha == -cone) then + do i=1, n_row + y(i) = -sv%d(i) * x(i) + end do + else + do i=1, n_row + y(i) = alpha * sv%d(i) * x(i) + end do + end if + + else if (beta == cone) then + + if (alpha == czero) then + !y(1:n_row) = czero + else if (alpha == cone) then + do i=1, n_row + y(i) = sv%d(i) * x(i) + y(i) + end do + else if (alpha == -cone) then + do i=1, n_row + y(i) = -sv%d(i) * x(i) + y(i) + end do + else + do i=1, n_row + y(i) = alpha * sv%d(i) * x(i) + y(i) + end do + end if + + else if (beta == -cone) then + + if (alpha == czero) then + y(1:n_row) = -y(1:n_row) + else if (alpha == cone) then + do i=1, n_row + y(i) = sv%d(i) * x(i) - y(i) + end do + else if (alpha == -cone) then + do i=1, n_row + y(i) = -sv%d(i) * x(i) - y(i) + end do + else + do i=1, n_row + y(i) = alpha * sv%d(i) * x(i) - y(i) + end do + end if + + else + + if (alpha == czero) then + y(1:n_row) = beta *y(1:n_row) + else if (alpha == cone) then + do i=1, n_row + y(i) = sv%d(i) * x(i) + beta*y(i) + end do + else if (alpha == -cone) then + do i=1, n_row + y(i) = -sv%d(i) * x(i) + beta*y(i) + end do + else + do i=1, n_row + y(i) = alpha * sv%d(i) * x(i) + beta*y(i) + end do + end if + + end if + + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + +end subroutine mld_c_diag_solver_apply + +subroutine mld_c_diag_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold) + use psb_base_mod + use mld_c_diag_solver, mld_protect_name => mld_c_diag_solver_bld + + Implicit None + + ! Arguments + type(psb_cspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_c_diag_solver_type), intent(inout) :: sv + character, intent(in) :: upd + integer, intent(out) :: info + type(psb_cspmat_type), intent(in), target, optional :: b + class(psb_c_base_sparse_mat), intent(in), optional :: amold + class(psb_c_base_vect_type), intent(in), optional :: vmold + ! Local variables + integer :: n_row,n_col, nrow_a, nztota + complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) + integer :: ictxt,np,me,i, err_act, debug_unit, debug_level + character(len=20) :: name='c_diag_solver_bld', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ictxt = desc_a%get_context() + call psb_info(ictxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + + n_row = desc_a%get_local_rows() + nrow_a = a%get_nrows() + if (allocated(sv%d)) then + if (size(sv%d) < n_row) then + deallocate(sv%d) + endif + endif + if (.not.allocated(sv%d)) then + allocate(sv%d(n_row),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + + endif + + call a%get_diag(sv%d,info) + if (present(b)) then + if (info == psb_success_) call b%get_diag(sv%d(nrow_a+1:), info) + end if + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='get_diag') + goto 9999 + end if + + do i=1,n_row + if (sv%d(i) == czero) then + sv%d(i) = cone + else + sv%d(i) = cone/sv%d(i) + end if + end do + allocate(sv%dv,stat=info) + if (info == psb_success_) then + if (present(vmold)) then + allocate(sv%dv%v,mold=vmold,stat=info) + else + allocate(psb_c_base_vect_type :: sv%dv%v,stat=info) + end if + end if + if (info == psb_success_) then + call sv%dv%bld(sv%d) + else + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate sv%dv') + goto 9999 + end if + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return +end subroutine mld_c_diag_solver_bld + diff --git a/mlprec/impl/mld_c_id_solver_impl.f90 b/mlprec/impl/mld_c_id_solver_impl.f90 new file mode 100644 index 00000000..4b24b3f9 --- /dev/null +++ b/mlprec/impl/mld_c_id_solver_impl.f90 @@ -0,0 +1,139 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.0) +!!$ +!!$ (C) Copyright 2008,2009,2010, 2010 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ 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 MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific 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 MLD2P4 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. +!!$ +!!$ +! +! +! +! Identity solver. Reference for nullprec. +! +! + + subroutine mld_c_id_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) + use psb_base_mod + use mld_c_id_solver, mld_protect_name => mld_c_id_solver_apply_vect + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_c_id_solver_type), intent(inout) :: sv + type(psb_c_vect_type),intent(inout) :: x + type(psb_c_vect_type),intent(inout) :: y + complex(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_spk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + + integer :: n_row,n_col + integer :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='c_id_solver_apply_vect' + + call psb_erractionsave(err_act) + + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T') + case('C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + call psb_geaxpby(alpha,x,beta,y,desc_data,info) + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine mld_c_id_solver_apply_vect + + + subroutine mld_c_id_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) + use psb_base_mod + use mld_c_id_solver, mld_protect_name => mld_c_id_solver_apply + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_c_id_solver_type), intent(in) :: sv + complex(psb_spk_),intent(inout) :: x(:) + complex(psb_spk_),intent(inout) :: y(:) + complex(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_spk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + + integer :: n_row,n_col + integer :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='c_id_solver_apply' + + call psb_erractionsave(err_act) + + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T') + case('C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + call psb_geaxpby(alpha,x,beta,y,desc_data,info) + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine mld_c_id_solver_apply diff --git a/mlprec/impl/mld_d_diag_solver_impl.f90 b/mlprec/impl/mld_d_diag_solver_impl.f90 new file mode 100644 index 00000000..0939865a --- /dev/null +++ b/mlprec/impl/mld_d_diag_solver_impl.f90 @@ -0,0 +1,413 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.0) +!!$ +!!$ (C) Copyright 2008,2009,2010, 2010 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ 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 MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific 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 MLD2P4 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 mld_d_diag_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) + use psb_base_mod + use mld_d_diag_solver, mld_protect_name => mld_d_diag_solver_apply_vect + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_d_diag_solver_type), intent(inout) :: sv + type(psb_d_vect_type), intent(inout) :: x + type(psb_d_vect_type), intent(inout) :: y + real(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + + integer :: n_row,n_col + real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) + integer :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='d_diag_solver_apply' + + call psb_erractionsave(err_act) + + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T','C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + n_row = desc_data%get_local_rows() + n_col = desc_data%get_local_cols() + if (x%get_nrows() < n_row) then + info = 36 + call psb_errpush(info,name,i_err=(/2,n_row,0,0,0/)) + goto 9999 + end if + if (y%get_nrows() < n_row) then + info = 36 + call psb_errpush(info,name,i_err=(/3,n_row,0,0,0/)) + goto 9999 + end if + if (.not.allocated(sv%dv)) then + info = 1124 + call psb_errpush(info,name,a_err="preconditioner: D") + goto 9999 + end if + if (sv%dv%get_nrows() < n_row) then + info = 1124 + call psb_errpush(info,name,a_err="preconditioner: D") + goto 9999 + end if + + call y%mlt(alpha,sv%dv,x,beta,info,conjgx=trans_) + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='vect%mlt') + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + +end subroutine mld_d_diag_solver_apply_vect + +subroutine mld_d_diag_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) + use psb_base_mod + use mld_d_diag_solver, mld_protect_name => mld_d_diag_solver_apply + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_d_diag_solver_type), intent(in) :: sv + real(psb_dpk_), intent(inout) :: x(:) + real(psb_dpk_), intent(inout) :: y(:) + real(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + + integer :: n_row,n_col + real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) + integer :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='d_diag_solver_apply' + + call psb_erractionsave(err_act) + + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T','C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + n_row = desc_data%get_local_rows() + n_col = desc_data%get_local_cols() + + if (trans_ == 'C') then + if (beta == dzero) then + + if (alpha == dzero) then + y(1:n_row) = dzero + else if (alpha == done) then + do i=1, n_row + y(i) = (sv%d(i)) * x(i) + end do + else if (alpha == -done) then + do i=1, n_row + y(i) = -(sv%d(i)) * x(i) + end do + else + do i=1, n_row + y(i) = alpha * (sv%d(i)) * x(i) + end do + end if + + else if (beta == done) then + + if (alpha == dzero) then + !y(1:n_row) = dzero + else if (alpha == done) then + do i=1, n_row + y(i) = (sv%d(i)) * x(i) + y(i) + end do + else if (alpha == -done) then + do i=1, n_row + y(i) = -(sv%d(i)) * x(i) + y(i) + end do + else + do i=1, n_row + y(i) = alpha * (sv%d(i)) * x(i) + y(i) + end do + end if + + else if (beta == -done) then + + if (alpha == dzero) then + y(1:n_row) = -y(1:n_row) + else if (alpha == done) then + do i=1, n_row + y(i) = (sv%d(i)) * x(i) - y(i) + end do + else if (alpha == -done) then + do i=1, n_row + y(i) = -(sv%d(i)) * x(i) - y(i) + end do + else + do i=1, n_row + y(i) = alpha * (sv%d(i)) * x(i) - y(i) + end do + end if + + else + + if (alpha == dzero) then + y(1:n_row) = beta *y(1:n_row) + else if (alpha == done) then + do i=1, n_row + y(i) = (sv%d(i)) * x(i) + beta*y(i) + end do + else if (alpha == -done) then + do i=1, n_row + y(i) = -(sv%d(i)) * x(i) + beta*y(i) + end do + else + do i=1, n_row + y(i) = alpha * (sv%d(i)) * x(i) + beta*y(i) + end do + end if + + end if + + else if (trans_ /= 'C') then + + if (beta == dzero) then + + if (alpha == dzero) then + y(1:n_row) = dzero + else if (alpha == done) then + do i=1, n_row + y(i) = sv%d(i) * x(i) + end do + else if (alpha == -done) then + do i=1, n_row + y(i) = -sv%d(i) * x(i) + end do + else + do i=1, n_row + y(i) = alpha * sv%d(i) * x(i) + end do + end if + + else if (beta == done) then + + if (alpha == dzero) then + !y(1:n_row) = dzero + else if (alpha == done) then + do i=1, n_row + y(i) = sv%d(i) * x(i) + y(i) + end do + else if (alpha == -done) then + do i=1, n_row + y(i) = -sv%d(i) * x(i) + y(i) + end do + else + do i=1, n_row + y(i) = alpha * sv%d(i) * x(i) + y(i) + end do + end if + + else if (beta == -done) then + + if (alpha == dzero) then + y(1:n_row) = -y(1:n_row) + else if (alpha == done) then + do i=1, n_row + y(i) = sv%d(i) * x(i) - y(i) + end do + else if (alpha == -done) then + do i=1, n_row + y(i) = -sv%d(i) * x(i) - y(i) + end do + else + do i=1, n_row + y(i) = alpha * sv%d(i) * x(i) - y(i) + end do + end if + + else + + if (alpha == dzero) then + y(1:n_row) = beta *y(1:n_row) + else if (alpha == done) then + do i=1, n_row + y(i) = sv%d(i) * x(i) + beta*y(i) + end do + else if (alpha == -done) then + do i=1, n_row + y(i) = -sv%d(i) * x(i) + beta*y(i) + end do + else + do i=1, n_row + y(i) = alpha * sv%d(i) * x(i) + beta*y(i) + end do + end if + + end if + + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + +end subroutine mld_d_diag_solver_apply + +subroutine mld_d_diag_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold) + use psb_base_mod + use mld_d_diag_solver, mld_protect_name => mld_d_diag_solver_bld + + Implicit None + + ! Arguments + type(psb_dspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_d_diag_solver_type), intent(inout) :: sv + character, intent(in) :: upd + integer, intent(out) :: info + type(psb_dspmat_type), intent(in), target, optional :: b + class(psb_d_base_sparse_mat), intent(in), optional :: amold + class(psb_d_base_vect_type), intent(in), optional :: vmold + ! Local variables + integer :: n_row,n_col, nrow_a, nztota + real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) + integer :: ictxt,np,me,i, err_act, debug_unit, debug_level + character(len=20) :: name='d_diag_solver_bld', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ictxt = desc_a%get_context() + call psb_info(ictxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + + n_row = desc_a%get_local_rows() + nrow_a = a%get_nrows() + if (allocated(sv%d)) then + if (size(sv%d) < n_row) then + deallocate(sv%d) + endif + endif + if (.not.allocated(sv%d)) then + allocate(sv%d(n_row),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + + endif + + call a%get_diag(sv%d,info) + if (present(b)) then + if (info == psb_success_) call b%get_diag(sv%d(nrow_a+1:), info) + end if + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='get_diag') + goto 9999 + end if + + do i=1,n_row + if (sv%d(i) == dzero) then + sv%d(i) = done + else + sv%d(i) = done/sv%d(i) + end if + end do + allocate(sv%dv,stat=info) + if (info == psb_success_) then + if (present(vmold)) then + allocate(sv%dv%v,mold=vmold,stat=info) + else + allocate(psb_d_base_vect_type :: sv%dv%v,stat=info) + end if + end if + if (info == psb_success_) then + call sv%dv%bld(sv%d) + else + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate sv%dv') + goto 9999 + end if + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return +end subroutine mld_d_diag_solver_bld + diff --git a/mlprec/impl/mld_d_id_solver_impl.f90 b/mlprec/impl/mld_d_id_solver_impl.f90 new file mode 100644 index 00000000..84c37136 --- /dev/null +++ b/mlprec/impl/mld_d_id_solver_impl.f90 @@ -0,0 +1,139 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.0) +!!$ +!!$ (C) Copyright 2008,2009,2010, 2010 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ 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 MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific 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 MLD2P4 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. +!!$ +!!$ +! +! +! +! Identity solver. Reference for nullprec. +! +! + + subroutine mld_d_id_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) + use psb_base_mod + use mld_d_id_solver, mld_protect_name => mld_d_id_solver_apply_vect + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_d_id_solver_type), intent(inout) :: sv + type(psb_d_vect_type),intent(inout) :: x + type(psb_d_vect_type),intent(inout) :: y + real(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + + integer :: n_row,n_col + integer :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='d_id_solver_apply_vect' + + call psb_erractionsave(err_act) + + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T') + case('C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + call psb_geaxpby(alpha,x,beta,y,desc_data,info) + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine mld_d_id_solver_apply_vect + + + subroutine mld_d_id_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) + use psb_base_mod + use mld_d_id_solver, mld_protect_name => mld_d_id_solver_apply + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_d_id_solver_type), intent(in) :: sv + real(psb_dpk_),intent(inout) :: x(:) + real(psb_dpk_),intent(inout) :: y(:) + real(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + + integer :: n_row,n_col + integer :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='d_id_solver_apply' + + call psb_erractionsave(err_act) + + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T') + case('C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + call psb_geaxpby(alpha,x,beta,y,desc_data,info) + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine mld_d_id_solver_apply diff --git a/mlprec/impl/mld_s_diag_solver_impl.f90 b/mlprec/impl/mld_s_diag_solver_impl.f90 new file mode 100644 index 00000000..42e7738e --- /dev/null +++ b/mlprec/impl/mld_s_diag_solver_impl.f90 @@ -0,0 +1,413 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.0) +!!$ +!!$ (C) Copyright 2008,2009,2010, 2010 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ 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 MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific 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 MLD2P4 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 mld_s_diag_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) + use psb_base_mod + use mld_s_diag_solver, mld_protect_name => mld_s_diag_solver_apply_vect + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_s_diag_solver_type), intent(inout) :: sv + type(psb_s_vect_type), intent(inout) :: x + type(psb_s_vect_type), intent(inout) :: y + real(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_spk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + + integer :: n_row,n_col + real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) + integer :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='s_diag_solver_apply' + + call psb_erractionsave(err_act) + + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T','C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + n_row = desc_data%get_local_rows() + n_col = desc_data%get_local_cols() + if (x%get_nrows() < n_row) then + info = 36 + call psb_errpush(info,name,i_err=(/2,n_row,0,0,0/)) + goto 9999 + end if + if (y%get_nrows() < n_row) then + info = 36 + call psb_errpush(info,name,i_err=(/3,n_row,0,0,0/)) + goto 9999 + end if + if (.not.allocated(sv%dv)) then + info = 1124 + call psb_errpush(info,name,a_err="preconditioner: D") + goto 9999 + end if + if (sv%dv%get_nrows() < n_row) then + info = 1124 + call psb_errpush(info,name,a_err="preconditioner: D") + goto 9999 + end if + + call y%mlt(alpha,sv%dv,x,beta,info,conjgx=trans_) + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='vect%mlt') + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + +end subroutine mld_s_diag_solver_apply_vect + +subroutine mld_s_diag_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) + use psb_base_mod + use mld_s_diag_solver, mld_protect_name => mld_s_diag_solver_apply + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_s_diag_solver_type), intent(in) :: sv + real(psb_spk_), intent(inout) :: x(:) + real(psb_spk_), intent(inout) :: y(:) + real(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_spk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + + integer :: n_row,n_col + real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) + integer :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='s_diag_solver_apply' + + call psb_erractionsave(err_act) + + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T','C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + n_row = desc_data%get_local_rows() + n_col = desc_data%get_local_cols() + + if (trans_ == 'C') then + if (beta == szero) then + + if (alpha == szero) then + y(1:n_row) = szero + else if (alpha == sone) then + do i=1, n_row + y(i) = (sv%d(i)) * x(i) + end do + else if (alpha == -sone) then + do i=1, n_row + y(i) = -(sv%d(i)) * x(i) + end do + else + do i=1, n_row + y(i) = alpha * (sv%d(i)) * x(i) + end do + end if + + else if (beta == sone) then + + if (alpha == szero) then + !y(1:n_row) = szero + else if (alpha == sone) then + do i=1, n_row + y(i) = (sv%d(i)) * x(i) + y(i) + end do + else if (alpha == -sone) then + do i=1, n_row + y(i) = -(sv%d(i)) * x(i) + y(i) + end do + else + do i=1, n_row + y(i) = alpha * (sv%d(i)) * x(i) + y(i) + end do + end if + + else if (beta == -sone) then + + if (alpha == szero) then + y(1:n_row) = -y(1:n_row) + else if (alpha == sone) then + do i=1, n_row + y(i) = (sv%d(i)) * x(i) - y(i) + end do + else if (alpha == -sone) then + do i=1, n_row + y(i) = -(sv%d(i)) * x(i) - y(i) + end do + else + do i=1, n_row + y(i) = alpha * (sv%d(i)) * x(i) - y(i) + end do + end if + + else + + if (alpha == szero) then + y(1:n_row) = beta *y(1:n_row) + else if (alpha == sone) then + do i=1, n_row + y(i) = (sv%d(i)) * x(i) + beta*y(i) + end do + else if (alpha == -sone) then + do i=1, n_row + y(i) = -(sv%d(i)) * x(i) + beta*y(i) + end do + else + do i=1, n_row + y(i) = alpha * (sv%d(i)) * x(i) + beta*y(i) + end do + end if + + end if + + else if (trans_ /= 'C') then + + if (beta == szero) then + + if (alpha == szero) then + y(1:n_row) = szero + else if (alpha == sone) then + do i=1, n_row + y(i) = sv%d(i) * x(i) + end do + else if (alpha == -sone) then + do i=1, n_row + y(i) = -sv%d(i) * x(i) + end do + else + do i=1, n_row + y(i) = alpha * sv%d(i) * x(i) + end do + end if + + else if (beta == sone) then + + if (alpha == szero) then + !y(1:n_row) = szero + else if (alpha == sone) then + do i=1, n_row + y(i) = sv%d(i) * x(i) + y(i) + end do + else if (alpha == -sone) then + do i=1, n_row + y(i) = -sv%d(i) * x(i) + y(i) + end do + else + do i=1, n_row + y(i) = alpha * sv%d(i) * x(i) + y(i) + end do + end if + + else if (beta == -sone) then + + if (alpha == szero) then + y(1:n_row) = -y(1:n_row) + else if (alpha == sone) then + do i=1, n_row + y(i) = sv%d(i) * x(i) - y(i) + end do + else if (alpha == -sone) then + do i=1, n_row + y(i) = -sv%d(i) * x(i) - y(i) + end do + else + do i=1, n_row + y(i) = alpha * sv%d(i) * x(i) - y(i) + end do + end if + + else + + if (alpha == szero) then + y(1:n_row) = beta *y(1:n_row) + else if (alpha == sone) then + do i=1, n_row + y(i) = sv%d(i) * x(i) + beta*y(i) + end do + else if (alpha == -sone) then + do i=1, n_row + y(i) = -sv%d(i) * x(i) + beta*y(i) + end do + else + do i=1, n_row + y(i) = alpha * sv%d(i) * x(i) + beta*y(i) + end do + end if + + end if + + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + +end subroutine mld_s_diag_solver_apply + +subroutine mld_s_diag_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold) + use psb_base_mod + use mld_s_diag_solver, mld_protect_name => mld_s_diag_solver_bld + + Implicit None + + ! Arguments + type(psb_sspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_s_diag_solver_type), intent(inout) :: sv + character, intent(in) :: upd + integer, intent(out) :: info + type(psb_sspmat_type), intent(in), target, optional :: b + class(psb_s_base_sparse_mat), intent(in), optional :: amold + class(psb_s_base_vect_type), intent(in), optional :: vmold + ! Local variables + integer :: n_row,n_col, nrow_a, nztota + real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) + integer :: ictxt,np,me,i, err_act, debug_unit, debug_level + character(len=20) :: name='s_diag_solver_bld', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ictxt = desc_a%get_context() + call psb_info(ictxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + + n_row = desc_a%get_local_rows() + nrow_a = a%get_nrows() + if (allocated(sv%d)) then + if (size(sv%d) < n_row) then + deallocate(sv%d) + endif + endif + if (.not.allocated(sv%d)) then + allocate(sv%d(n_row),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + + endif + + call a%get_diag(sv%d,info) + if (present(b)) then + if (info == psb_success_) call b%get_diag(sv%d(nrow_a+1:), info) + end if + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='get_diag') + goto 9999 + end if + + do i=1,n_row + if (sv%d(i) == szero) then + sv%d(i) = sone + else + sv%d(i) = sone/sv%d(i) + end if + end do + allocate(sv%dv,stat=info) + if (info == psb_success_) then + if (present(vmold)) then + allocate(sv%dv%v,mold=vmold,stat=info) + else + allocate(psb_s_base_vect_type :: sv%dv%v,stat=info) + end if + end if + if (info == psb_success_) then + call sv%dv%bld(sv%d) + else + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate sv%dv') + goto 9999 + end if + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return +end subroutine mld_s_diag_solver_bld + diff --git a/mlprec/impl/mld_s_id_solver_impl.f90 b/mlprec/impl/mld_s_id_solver_impl.f90 new file mode 100644 index 00000000..2420f491 --- /dev/null +++ b/mlprec/impl/mld_s_id_solver_impl.f90 @@ -0,0 +1,139 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.0) +!!$ +!!$ (C) Copyright 2008,2009,2010, 2010 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ 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 MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific 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 MLD2P4 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. +!!$ +!!$ +! +! +! +! Identity solver. Reference for nullprec. +! +! + + subroutine mld_s_id_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) + use psb_base_mod + use mld_s_id_solver, mld_protect_name => mld_s_id_solver_apply_vect + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_s_id_solver_type), intent(inout) :: sv + type(psb_s_vect_type),intent(inout) :: x + type(psb_s_vect_type),intent(inout) :: y + real(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_spk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + + integer :: n_row,n_col + integer :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='s_id_solver_apply_vect' + + call psb_erractionsave(err_act) + + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T') + case('C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + call psb_geaxpby(alpha,x,beta,y,desc_data,info) + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine mld_s_id_solver_apply_vect + + + subroutine mld_s_id_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) + use psb_base_mod + use mld_s_id_solver, mld_protect_name => mld_s_id_solver_apply + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_s_id_solver_type), intent(in) :: sv + real(psb_spk_),intent(inout) :: x(:) + real(psb_spk_),intent(inout) :: y(:) + real(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_spk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + + integer :: n_row,n_col + integer :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='s_id_solver_apply' + + call psb_erractionsave(err_act) + + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T') + case('C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + call psb_geaxpby(alpha,x,beta,y,desc_data,info) + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine mld_s_id_solver_apply diff --git a/mlprec/impl/mld_z_diag_solver_impl.f90 b/mlprec/impl/mld_z_diag_solver_impl.f90 new file mode 100644 index 00000000..190da178 --- /dev/null +++ b/mlprec/impl/mld_z_diag_solver_impl.f90 @@ -0,0 +1,413 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.0) +!!$ +!!$ (C) Copyright 2008,2009,2010, 2010 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ 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 MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific 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 MLD2P4 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 mld_z_diag_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) + use psb_base_mod + use mld_z_diag_solver, mld_protect_name => mld_z_diag_solver_apply_vect + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_z_diag_solver_type), intent(inout) :: sv + type(psb_z_vect_type), intent(inout) :: x + type(psb_z_vect_type), intent(inout) :: y + complex(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + + integer :: n_row,n_col + complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) + integer :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='z_diag_solver_apply' + + call psb_erractionsave(err_act) + + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T','C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + n_row = desc_data%get_local_rows() + n_col = desc_data%get_local_cols() + if (x%get_nrows() < n_row) then + info = 36 + call psb_errpush(info,name,i_err=(/2,n_row,0,0,0/)) + goto 9999 + end if + if (y%get_nrows() < n_row) then + info = 36 + call psb_errpush(info,name,i_err=(/3,n_row,0,0,0/)) + goto 9999 + end if + if (.not.allocated(sv%dv)) then + info = 1124 + call psb_errpush(info,name,a_err="preconditioner: D") + goto 9999 + end if + if (sv%dv%get_nrows() < n_row) then + info = 1124 + call psb_errpush(info,name,a_err="preconditioner: D") + goto 9999 + end if + + call y%mlt(alpha,sv%dv,x,beta,info,conjgx=trans_) + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='vect%mlt') + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + +end subroutine mld_z_diag_solver_apply_vect + +subroutine mld_z_diag_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) + use psb_base_mod + use mld_z_diag_solver, mld_protect_name => mld_z_diag_solver_apply + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_z_diag_solver_type), intent(in) :: sv + complex(psb_dpk_), intent(inout) :: x(:) + complex(psb_dpk_), intent(inout) :: y(:) + complex(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + + integer :: n_row,n_col + complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) + integer :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='z_diag_solver_apply' + + call psb_erractionsave(err_act) + + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T','C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + n_row = desc_data%get_local_rows() + n_col = desc_data%get_local_cols() + + if (trans_ == 'C') then + if (beta == zzero) then + + if (alpha == zzero) then + y(1:n_row) = zzero + else if (alpha == zone) then + do i=1, n_row + y(i) = conjg(sv%d(i)) * x(i) + end do + else if (alpha == -zone) then + do i=1, n_row + y(i) = -conjg(sv%d(i)) * x(i) + end do + else + do i=1, n_row + y(i) = alpha * conjg(sv%d(i)) * x(i) + end do + end if + + else if (beta == zone) then + + if (alpha == zzero) then + !y(1:n_row) = zzero + else if (alpha == zone) then + do i=1, n_row + y(i) = conjg(sv%d(i)) * x(i) + y(i) + end do + else if (alpha == -zone) then + do i=1, n_row + y(i) = -conjg(sv%d(i)) * x(i) + y(i) + end do + else + do i=1, n_row + y(i) = alpha * conjg(sv%d(i)) * x(i) + y(i) + end do + end if + + else if (beta == -zone) then + + if (alpha == zzero) then + y(1:n_row) = -y(1:n_row) + else if (alpha == zone) then + do i=1, n_row + y(i) = conjg(sv%d(i)) * x(i) - y(i) + end do + else if (alpha == -zone) then + do i=1, n_row + y(i) = -conjg(sv%d(i)) * x(i) - y(i) + end do + else + do i=1, n_row + y(i) = alpha * conjg(sv%d(i)) * x(i) - y(i) + end do + end if + + else + + if (alpha == zzero) then + y(1:n_row) = beta *y(1:n_row) + else if (alpha == zone) then + do i=1, n_row + y(i) = conjg(sv%d(i)) * x(i) + beta*y(i) + end do + else if (alpha == -zone) then + do i=1, n_row + y(i) = -conjg(sv%d(i)) * x(i) + beta*y(i) + end do + else + do i=1, n_row + y(i) = alpha * conjg(sv%d(i)) * x(i) + beta*y(i) + end do + end if + + end if + + else if (trans_ /= 'C') then + + if (beta == zzero) then + + if (alpha == zzero) then + y(1:n_row) = zzero + else if (alpha == zone) then + do i=1, n_row + y(i) = sv%d(i) * x(i) + end do + else if (alpha == -zone) then + do i=1, n_row + y(i) = -sv%d(i) * x(i) + end do + else + do i=1, n_row + y(i) = alpha * sv%d(i) * x(i) + end do + end if + + else if (beta == zone) then + + if (alpha == zzero) then + !y(1:n_row) = zzero + else if (alpha == zone) then + do i=1, n_row + y(i) = sv%d(i) * x(i) + y(i) + end do + else if (alpha == -zone) then + do i=1, n_row + y(i) = -sv%d(i) * x(i) + y(i) + end do + else + do i=1, n_row + y(i) = alpha * sv%d(i) * x(i) + y(i) + end do + end if + + else if (beta == -zone) then + + if (alpha == zzero) then + y(1:n_row) = -y(1:n_row) + else if (alpha == zone) then + do i=1, n_row + y(i) = sv%d(i) * x(i) - y(i) + end do + else if (alpha == -zone) then + do i=1, n_row + y(i) = -sv%d(i) * x(i) - y(i) + end do + else + do i=1, n_row + y(i) = alpha * sv%d(i) * x(i) - y(i) + end do + end if + + else + + if (alpha == zzero) then + y(1:n_row) = beta *y(1:n_row) + else if (alpha == zone) then + do i=1, n_row + y(i) = sv%d(i) * x(i) + beta*y(i) + end do + else if (alpha == -zone) then + do i=1, n_row + y(i) = -sv%d(i) * x(i) + beta*y(i) + end do + else + do i=1, n_row + y(i) = alpha * sv%d(i) * x(i) + beta*y(i) + end do + end if + + end if + + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + +end subroutine mld_z_diag_solver_apply + +subroutine mld_z_diag_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold) + use psb_base_mod + use mld_z_diag_solver, mld_protect_name => mld_z_diag_solver_bld + + Implicit None + + ! Arguments + type(psb_zspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_z_diag_solver_type), intent(inout) :: sv + character, intent(in) :: upd + integer, intent(out) :: info + type(psb_zspmat_type), intent(in), target, optional :: b + class(psb_z_base_sparse_mat), intent(in), optional :: amold + class(psb_z_base_vect_type), intent(in), optional :: vmold + ! Local variables + integer :: n_row,n_col, nrow_a, nztota + complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) + integer :: ictxt,np,me,i, err_act, debug_unit, debug_level + character(len=20) :: name='z_diag_solver_bld', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ictxt = desc_a%get_context() + call psb_info(ictxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + + n_row = desc_a%get_local_rows() + nrow_a = a%get_nrows() + if (allocated(sv%d)) then + if (size(sv%d) < n_row) then + deallocate(sv%d) + endif + endif + if (.not.allocated(sv%d)) then + allocate(sv%d(n_row),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + end if + + endif + + call a%get_diag(sv%d,info) + if (present(b)) then + if (info == psb_success_) call b%get_diag(sv%d(nrow_a+1:), info) + end if + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='get_diag') + goto 9999 + end if + + do i=1,n_row + if (sv%d(i) == zzero) then + sv%d(i) = zone + else + sv%d(i) = zone/sv%d(i) + end if + end do + allocate(sv%dv,stat=info) + if (info == psb_success_) then + if (present(vmold)) then + allocate(sv%dv%v,mold=vmold,stat=info) + else + allocate(psb_z_base_vect_type :: sv%dv%v,stat=info) + end if + end if + if (info == psb_success_) then + call sv%dv%bld(sv%d) + else + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate sv%dv') + goto 9999 + end if + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return +end subroutine mld_z_diag_solver_bld + diff --git a/mlprec/impl/mld_z_id_solver_impl.f90 b/mlprec/impl/mld_z_id_solver_impl.f90 new file mode 100644 index 00000000..fee79225 --- /dev/null +++ b/mlprec/impl/mld_z_id_solver_impl.f90 @@ -0,0 +1,139 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.0) +!!$ +!!$ (C) Copyright 2008,2009,2010, 2010 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ 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 MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific 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 MLD2P4 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. +!!$ +!!$ +! +! +! +! Identity solver. Reference for nullprec. +! +! + + subroutine mld_z_id_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) + use psb_base_mod + use mld_z_id_solver, mld_protect_name => mld_z_id_solver_apply_vect + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_z_id_solver_type), intent(inout) :: sv + type(psb_z_vect_type),intent(inout) :: x + type(psb_z_vect_type),intent(inout) :: y + complex(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + + integer :: n_row,n_col + integer :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='z_id_solver_apply_vect' + + call psb_erractionsave(err_act) + + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T') + case('C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + call psb_geaxpby(alpha,x,beta,y,desc_data,info) + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine mld_z_id_solver_apply_vect + + + subroutine mld_z_id_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) + use psb_base_mod + use mld_z_id_solver, mld_protect_name => mld_z_id_solver_apply + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_z_id_solver_type), intent(in) :: sv + complex(psb_dpk_),intent(inout) :: x(:) + complex(psb_dpk_),intent(inout) :: y(:) + complex(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + + integer :: n_row,n_col + integer :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='z_id_solver_apply' + + call psb_erractionsave(err_act) + + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T') + case('C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + call psb_geaxpby(alpha,x,beta,y,desc_data,info) + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine mld_z_id_solver_apply diff --git a/mlprec/mld_c_diag_solver.f90 b/mlprec/mld_c_diag_solver.f90 index 8d68b864..cae39949 100644 --- a/mlprec/mld_c_diag_solver.f90 +++ b/mlprec/mld_c_diag_solver.f90 @@ -51,9 +51,9 @@ module mld_c_diag_solver type(psb_c_vect_type), allocatable :: dv complex(psb_spk_), allocatable :: d(:) contains - procedure, pass(sv) :: build => c_diag_solver_bld - procedure, pass(sv) :: apply_v => c_diag_solver_apply_vect - procedure, pass(sv) :: apply_a => c_diag_solver_apply + procedure, pass(sv) :: build => mld_c_diag_solver_bld + procedure, pass(sv) :: apply_v => mld_c_diag_solver_apply_vect + procedure, pass(sv) :: apply_a => mld_c_diag_solver_apply procedure, pass(sv) :: free => c_diag_solver_free procedure, pass(sv) :: seti => c_diag_solver_seti procedure, pass(sv) :: setc => c_diag_solver_setc @@ -64,385 +64,63 @@ module mld_c_diag_solver end type mld_c_diag_solver_type - private :: c_diag_solver_bld, c_diag_solver_apply, & - & c_diag_solver_free, c_diag_solver_seti, & + private :: c_diag_solver_free, c_diag_solver_seti, & & c_diag_solver_setc, c_diag_solver_setr,& & c_diag_solver_descr, c_diag_solver_sizeof,& & c_diag_solver_get_nzeros + interface mld_c_diag_solver_apply_vect + subroutine mld_c_diag_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) + import :: psb_desc_type, psb_cspmat_type, psb_c_base_sparse_mat, & + & psb_c_vect_type, psb_c_base_vect_type, psb_spk_, mld_c_diag_solver_type + type(psb_desc_type), intent(in) :: desc_data + class(mld_c_diag_solver_type), intent(inout) :: sv + type(psb_c_vect_type), intent(inout) :: x + type(psb_c_vect_type), intent(inout) :: y + complex(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_spk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + end subroutine mld_c_diag_solver_apply_vect + end interface mld_c_diag_solver_apply_vect + + interface mld_c_diag_solver_apply + subroutine mld_c_diag_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) + import :: psb_desc_type, psb_cspmat_type, psb_c_base_sparse_mat, & + & psb_c_vect_type, psb_c_base_vect_type, psb_spk_, mld_c_diag_solver_type + type(psb_desc_type), intent(in) :: desc_data + class(mld_c_diag_solver_type), intent(in) :: sv + complex(psb_spk_), intent(inout) :: x(:) + complex(psb_spk_), intent(inout) :: y(:) + complex(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_spk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + end subroutine mld_c_diag_solver_apply + end interface mld_c_diag_solver_apply + + interface mld_c_diag_solver_bld + subroutine mld_c_diag_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold) + import :: psb_desc_type, psb_cspmat_type, psb_c_base_sparse_mat, & + & psb_c_vect_type, psb_c_base_vect_type, psb_spk_, mld_c_diag_solver_type + type(psb_cspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_c_diag_solver_type), intent(inout) :: sv + character, intent(in) :: upd + integer, intent(out) :: info + type(psb_cspmat_type), intent(in), target, optional :: b + class(psb_c_base_sparse_mat), intent(in), optional :: amold + class(psb_c_base_vect_type), intent(in), optional :: vmold + end subroutine mld_c_diag_solver_bld + end interface mld_c_diag_solver_bld + + contains - subroutine c_diag_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) - use psb_base_mod - type(psb_desc_type), intent(in) :: desc_data - class(mld_c_diag_solver_type), intent(inout) :: sv - type(psb_c_vect_type), intent(inout) :: x - type(psb_c_vect_type), intent(inout) :: y - complex(psb_spk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - complex(psb_spk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - - integer :: n_row,n_col - complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act - character :: trans_ - character(len=20) :: name='c_diag_solver_apply' - - call psb_erractionsave(err_act) - - info = psb_success_ - - trans_ = psb_toupper(trans) - select case(trans_) - case('N') - case('T','C') - case default - call psb_errpush(psb_err_iarg_invalid_i_,name) - goto 9999 - end select - - n_row = desc_data%get_local_rows() - n_col = desc_data%get_local_cols() - if (x%get_nrows() < n_row) then - info = 36 - call psb_errpush(info,name,i_err=(/2,nrow,0,0,0/)) - goto 9999 - end if - if (y%get_nrows() < n_row) then - info = 36 - call psb_errpush(info,name,i_err=(/3,nrow,0,0,0/)) - goto 9999 - end if - if (.not.allocated(sv%dv)) then - info = 1124 - call psb_errpush(info,name,a_err="preconditioner: D") - goto 9999 - end if - if (sv%dv%get_nrows() < n_row) then - info = 1124 - call psb_errpush(info,name,a_err="preconditioner: D") - goto 9999 - end if - - call y%mlt(alpha,sv%dv,x,beta,info,trans=trans_) - - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='vect%mlt') - goto 9999 - end if - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - - end subroutine c_diag_solver_apply_vect - - subroutine c_diag_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) - use psb_base_mod - type(psb_desc_type), intent(in) :: desc_data - class(mld_c_diag_solver_type), intent(in) :: sv - complex(psb_spk_), intent(inout) :: x(:) - complex(psb_spk_), intent(inout) :: y(:) - complex(psb_spk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - complex(psb_spk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - - integer :: n_row,n_col - complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act - character :: trans_ - character(len=20) :: name='c_diag_solver_apply' - - call psb_erractionsave(err_act) - - info = psb_success_ - - trans_ = psb_toupper(trans) - select case(trans_) - case('N') - case('T','C') - case default - call psb_errpush(psb_err_iarg_invalid_i_,name) - goto 9999 - end select - - n_row = desc_data%get_local_rows() - n_col = desc_data%get_local_cols() - - if (trans_ == 'C') then - if (beta == czero) then - - if (alpha == czero) then - y(1:n_row) = czero - else if (alpha == cone) then - do i=1, n_row - y(i) = conjg(sv%d(i)) * x(i) - end do - else if (alpha == -cone) then - do i=1, n_row - y(i) = -conjg(sv%d(i)) * x(i) - end do - else - do i=1, n_row - y(i) = alpha * conjg(sv%d(i)) * x(i) - end do - end if - - else if (beta == cone) then - - if (alpha == czero) then - !y(1:n_row) = czero - else if (alpha == cone) then - do i=1, n_row - y(i) = conjg(sv%d(i)) * x(i) + y(i) - end do - else if (alpha == -cone) then - do i=1, n_row - y(i) = -conjg(sv%d(i)) * x(i) + y(i) - end do - else - do i=1, n_row - y(i) = alpha * conjg(sv%d(i)) * x(i) + y(i) - end do - end if - - else if (beta == -cone) then - - if (alpha == czero) then - y(1:n_row) = -y(1:n_row) - else if (alpha == cone) then - do i=1, n_row - y(i) = conjg(sv%d(i)) * x(i) - y(i) - end do - else if (alpha == -cone) then - do i=1, n_row - y(i) = -conjg(sv%d(i)) * x(i) - y(i) - end do - else - do i=1, n_row - y(i) = alpha * conjg(sv%d(i)) * x(i) - y(i) - end do - end if - - else - - if (alpha == czero) then - y(1:n_row) = beta *y(1:n_row) - else if (alpha == cone) then - do i=1, n_row - y(i) = conjg(sv%d(i)) * x(i) + beta*y(i) - end do - else if (alpha == -cone) then - do i=1, n_row - y(i) = -conjg(sv%d(i)) * x(i) + beta*y(i) - end do - else - do i=1, n_row - y(i) = alpha * conjg(sv%d(i)) * x(i) + beta*y(i) - end do - end if - - end if - - else if (trans_ /= 'C') then - - if (beta == czero) then - - if (alpha == czero) then - y(1:n_row) = czero - else if (alpha == cone) then - do i=1, n_row - y(i) = sv%d(i) * x(i) - end do - else if (alpha == -cone) then - do i=1, n_row - y(i) = -sv%d(i) * x(i) - end do - else - do i=1, n_row - y(i) = alpha * sv%d(i) * x(i) - end do - end if - - else if (beta == cone) then - - if (alpha == czero) then - !y(1:n_row) = czero - else if (alpha == cone) then - do i=1, n_row - y(i) = sv%d(i) * x(i) + y(i) - end do - else if (alpha == -cone) then - do i=1, n_row - y(i) = -sv%d(i) * x(i) + y(i) - end do - else - do i=1, n_row - y(i) = alpha * sv%d(i) * x(i) + y(i) - end do - end if - - else if (beta == -cone) then - - if (alpha == czero) then - y(1:n_row) = -y(1:n_row) - else if (alpha == cone) then - do i=1, n_row - y(i) = sv%d(i) * x(i) - y(i) - end do - else if (alpha == -cone) then - do i=1, n_row - y(i) = -sv%d(i) * x(i) - y(i) - end do - else - do i=1, n_row - y(i) = alpha * sv%d(i) * x(i) - y(i) - end do - end if - - else - - if (alpha == czero) then - y(1:n_row) = beta *y(1:n_row) - else if (alpha == cone) then - do i=1, n_row - y(i) = sv%d(i) * x(i) + beta*y(i) - end do - else if (alpha == -cone) then - do i=1, n_row - y(i) = -sv%d(i) * x(i) + beta*y(i) - end do - else - do i=1, n_row - y(i) = alpha * sv%d(i) * x(i) + beta*y(i) - end do - end if - - end if - - end if - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - - end subroutine c_diag_solver_apply - - subroutine c_diag_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold) - - use psb_base_mod - - Implicit None - - ! Arguments - type(psb_cspmat_type), intent(in), target :: a - Type(psb_desc_type), Intent(in) :: desc_a - class(mld_c_diag_solver_type), intent(inout) :: sv - character, intent(in) :: upd - integer, intent(out) :: info - type(psb_cspmat_type), intent(in), target, optional :: b - class(psb_c_base_sparse_mat), intent(in), optional :: amold - class(psb_c_base_vect_type), intent(in), optional :: vmold - ! Local variables - integer :: n_row,n_col, nrow_a, nztota - complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act, debug_unit, debug_level - character(len=20) :: name='c_diag_solver_bld', ch_err - - info=psb_success_ - call psb_erractionsave(err_act) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - ictxt = desc_a%get_context() - call psb_info(ictxt, me, np) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),' start' - - - n_row = desc_a%get_local_rows() - nrow_a = a%get_nrows() - if (allocated(sv%d)) then - if (size(sv%d) < n_row) then - deallocate(sv%d) - endif - endif - if (.not.allocated(sv%d)) then - allocate(sv%d(n_row),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 - end if - - endif - - call a%get_diag(sv%d,info) - if (present(b)) then - if (info == psb_success_) call b%get_diag(sv%d(nrow_a+1:), info) - end if - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='get_diag') - goto 9999 - end if - - do i=1,n_row - if (sv%d(i) == czero) then - sv%d(i) = cone - else - sv%d(i) = cone/sv%d(i) - end if - end do - allocate(sv%dv,stat=info) - if (info == psb_success_) then - if (present(vmold)) then - allocate(sv%dv%v,mold=vmold,stat=info) - else - allocate(psb_c_base_vect_type :: sv%dv%v,stat=info) - end if - end if - if (info == psb_success_) then - call sv%dv%bld(sv%d) - else - call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate sv%dv') - goto 9999 - end if - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),' end' - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - end subroutine c_diag_solver_bld - subroutine c_diag_solver_seti(sv,what,val,info) - use psb_base_mod - Implicit None ! Arguments @@ -454,34 +132,14 @@ contains character(len=20) :: name='c_diag_solver_seti' info = psb_success_ -!!$ call psb_erractionsave(err_act) -!!$ -!!$ select case(what) -!!$ case(mld_sub_solve_) -!!$ sv%fact_type = val -!!$ case(mld_sub_fillin_) -!!$ sv%fill_in = val -!!$ case default -!!$ write(0,*) name,': Error: invalid WHAT' -!!$ info = -2 -!!$ end select -!!$ -!!$ call psb_erractionrestore(err_act) - return -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if + return + end subroutine c_diag_solver_seti subroutine c_diag_solver_setc(sv,what,val,info) - use psb_base_mod - Implicit None ! Arguments @@ -493,33 +151,12 @@ contains character(len=20) :: name='c_diag_solver_setc' info = psb_success_ -!!$ call psb_erractionsave(err_act) -!!$ -!!$ -!!$ call mld_stringval(val,ival,info) -!!$ if (info == psb_success_) call sv%set(what,ival,info) -!!$ if (info /= psb_success_) then -!!$ info = psb_err_from_subroutine_ -!!$ call psb_errpush(info, name) -!!$ goto 9999 -!!$ end if -!!$ -!!$ call psb_erractionrestore(err_act) - return -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if return end subroutine c_diag_solver_setc subroutine c_diag_solver_setr(sv,what,val,info) - use psb_base_mod - Implicit None ! Arguments @@ -531,33 +168,13 @@ contains character(len=20) :: name='c_diag_solver_setr' info = psb_success_ -!!$ call psb_erractionsave(err_act) -!!$ -!!$ select case(what) -!!$ case(mld_sub_iluthrs_) -!!$ sv%thresh = val -!!$ case default -!!$ write(0,*) name,': Error: invalid WHAT' -!!$ info = -2 -!!$ goto 9999 -!!$ end select -!!$ -!!$ call psb_erractionrestore(err_act) - return -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if return + end subroutine c_diag_solver_setr subroutine c_diag_solver_free(sv,info) - use psb_base_mod - Implicit None ! Arguments @@ -594,8 +211,6 @@ contains subroutine c_diag_solver_descr(sv,info,iout,coarse) - use psb_base_mod - Implicit None ! Arguments @@ -624,7 +239,6 @@ contains end subroutine c_diag_solver_descr function c_diag_solver_sizeof(sv) result(val) - use psb_base_mod implicit none ! Arguments class(mld_c_diag_solver_type), intent(in) :: sv @@ -638,7 +252,6 @@ contains end function c_diag_solver_sizeof function c_diag_solver_get_nzeros(sv) result(val) - use psb_base_mod implicit none ! Arguments class(mld_c_diag_solver_type), intent(in) :: sv diff --git a/mlprec/mld_c_id_solver.f90 b/mlprec/mld_c_id_solver.f90 index 7fc411f4..5e5e654b 100644 --- a/mlprec/mld_c_id_solver.f90 +++ b/mlprec/mld_c_id_solver.f90 @@ -50,8 +50,8 @@ module mld_c_id_solver type, extends(mld_c_base_solver_type) :: mld_c_id_solver_type contains procedure, pass(sv) :: build => c_id_solver_bld - procedure, pass(sv) :: apply_v => c_id_solver_apply_vect - procedure, pass(sv) :: apply_a => c_id_solver_apply + procedure, pass(sv) :: apply_v => mld_c_id_solver_apply_vect + procedure, pass(sv) :: apply_a => mld_c_id_solver_apply procedure, pass(sv) :: free => c_id_solver_free procedure, pass(sv) :: seti => c_id_solver_seti procedure, pass(sv) :: setc => c_id_solver_setc @@ -60,111 +60,46 @@ module mld_c_id_solver end type mld_c_id_solver_type - private :: c_id_solver_bld, c_id_solver_apply, & + private :: c_id_solver_bld, & & c_id_solver_free, c_id_solver_seti, & & c_id_solver_setc, c_id_solver_setr,& & c_id_solver_descr + interface mld_c_id_solver_apply_vect + subroutine mld_c_id_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) + import :: psb_desc_type, psb_cspmat_type, psb_c_base_sparse_mat, & + & psb_c_vect_type, psb_c_base_vect_type, psb_spk_, mld_c_id_solver_type + type(psb_desc_type), intent(in) :: desc_data + class(mld_c_id_solver_type), intent(inout) :: sv + type(psb_c_vect_type),intent(inout) :: x + type(psb_c_vect_type),intent(inout) :: y + complex(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_spk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + end subroutine mld_c_id_solver_apply_vect + end interface mld_c_id_solver_apply_vect + + interface mld_c_id_solver_apply + subroutine mld_c_id_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) + import :: psb_desc_type, psb_cspmat_type, psb_c_base_sparse_mat, & + & psb_c_vect_type, psb_c_base_vect_type, psb_spk_, mld_c_id_solver_type + type(psb_desc_type), intent(in) :: desc_data + class(mld_c_id_solver_type), intent(in) :: sv + complex(psb_spk_),intent(inout) :: x(:) + complex(psb_spk_),intent(inout) :: y(:) + complex(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_spk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + end subroutine mld_c_id_solver_apply + end interface mld_c_id_solver_apply contains - subroutine c_id_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) - use psb_base_mod - type(psb_desc_type), intent(in) :: desc_data - class(mld_c_id_solver_type), intent(inout) :: sv - type(psb_c_vect_type),intent(inout) :: x - type(psb_c_vect_type),intent(inout) :: y - complex(psb_spk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - complex(psb_spk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - - integer :: n_row,n_col - complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act - character :: trans_ - character(len=20) :: name='c_id_solver_apply' - - call psb_erractionsave(err_act) - - info = psb_success_ - - trans_ = psb_toupper(trans) - select case(trans_) - case('N') - case('T') - case('C') - case default - call psb_errpush(psb_err_iarg_invalid_i_,name) - goto 9999 - end select - - call psb_geaxpby(alpha,x,beta,y,desc_data,info) - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - - end subroutine c_id_solver_apply_vect - - - subroutine c_id_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) - use psb_base_mod - type(psb_desc_type), intent(in) :: desc_data - class(mld_c_id_solver_type), intent(in) :: sv - complex(psb_spk_),intent(inout) :: x(:) - complex(psb_spk_),intent(inout) :: y(:) - complex(psb_spk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - complex(psb_spk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - - integer :: n_row,n_col - complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act - character :: trans_ - character(len=20) :: name='c_id_solver_apply' - - call psb_erractionsave(err_act) - - info = psb_success_ - - trans_ = psb_toupper(trans) - select case(trans_) - case('N') - case('T') - case('C') - case default - call psb_errpush(psb_err_iarg_invalid_i_,name) - goto 9999 - end select - - call psb_geaxpby(alpha,x,beta,y,desc_data,info) - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - - end subroutine c_id_solver_apply subroutine c_id_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold) - use psb_base_mod - Implicit None ! Arguments @@ -183,35 +118,13 @@ contains character(len=20) :: name='c_id_solver_bld', ch_err info=psb_success_ - call psb_erractionsave(err_act) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - ictxt = desc_a%get_context() - call psb_info(ictxt, me, np) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),' start' - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),' end' - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if return end subroutine c_id_solver_bld subroutine c_id_solver_seti(sv,what,val,info) - use psb_base_mod - Implicit None ! Arguments @@ -230,8 +143,6 @@ contains subroutine c_id_solver_setc(sv,what,val,info) - use psb_base_mod - Implicit None ! Arguments @@ -249,8 +160,6 @@ contains subroutine c_id_solver_setr(sv,what,val,info) - use psb_base_mod - Implicit None ! Arguments @@ -269,8 +178,6 @@ contains subroutine c_id_solver_free(sv,info) - use psb_base_mod - Implicit None ! Arguments @@ -286,8 +193,6 @@ contains subroutine c_id_solver_descr(sv,info,iout,coarse) - use psb_base_mod - Implicit None ! Arguments diff --git a/mlprec/mld_d_diag_solver.f90 b/mlprec/mld_d_diag_solver.f90 index 1cec1225..50f4c25e 100644 --- a/mlprec/mld_d_diag_solver.f90 +++ b/mlprec/mld_d_diag_solver.f90 @@ -51,9 +51,9 @@ module mld_d_diag_solver type(psb_d_vect_type), allocatable :: dv real(psb_dpk_), allocatable :: d(:) contains - procedure, pass(sv) :: build => d_diag_solver_bld - procedure, pass(sv) :: apply_v => d_diag_solver_apply_vect - procedure, pass(sv) :: apply_a => d_diag_solver_apply + procedure, pass(sv) :: build => mld_d_diag_solver_bld + procedure, pass(sv) :: apply_v => mld_d_diag_solver_apply_vect + procedure, pass(sv) :: apply_a => mld_d_diag_solver_apply procedure, pass(sv) :: free => d_diag_solver_free procedure, pass(sv) :: seti => d_diag_solver_seti procedure, pass(sv) :: setc => d_diag_solver_setc @@ -64,306 +64,63 @@ module mld_d_diag_solver end type mld_d_diag_solver_type - private :: d_diag_solver_bld, d_diag_solver_apply, & - & d_diag_solver_free, d_diag_solver_seti, & + private :: d_diag_solver_free, d_diag_solver_seti, & & d_diag_solver_setc, d_diag_solver_setr,& & d_diag_solver_descr, d_diag_solver_sizeof,& & d_diag_solver_get_nzeros + interface mld_d_diag_solver_apply_vect + subroutine mld_d_diag_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) + import :: psb_desc_type, psb_dspmat_type, psb_d_base_sparse_mat, & + & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, mld_d_diag_solver_type + type(psb_desc_type), intent(in) :: desc_data + class(mld_d_diag_solver_type), intent(inout) :: sv + type(psb_d_vect_type), intent(inout) :: x + type(psb_d_vect_type), intent(inout) :: y + real(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + end subroutine mld_d_diag_solver_apply_vect + end interface mld_d_diag_solver_apply_vect + + interface mld_d_diag_solver_apply + subroutine mld_d_diag_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) + import :: psb_desc_type, psb_dspmat_type, psb_d_base_sparse_mat, & + & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, mld_d_diag_solver_type + type(psb_desc_type), intent(in) :: desc_data + class(mld_d_diag_solver_type), intent(in) :: sv + real(psb_dpk_), intent(inout) :: x(:) + real(psb_dpk_), intent(inout) :: y(:) + real(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + end subroutine mld_d_diag_solver_apply + end interface mld_d_diag_solver_apply + + interface mld_d_diag_solver_bld + subroutine mld_d_diag_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold) + import :: psb_desc_type, psb_dspmat_type, psb_d_base_sparse_mat, & + & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, mld_d_diag_solver_type + type(psb_dspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_d_diag_solver_type), intent(inout) :: sv + character, intent(in) :: upd + integer, intent(out) :: info + type(psb_dspmat_type), intent(in), target, optional :: b + class(psb_d_base_sparse_mat), intent(in), optional :: amold + class(psb_d_base_vect_type), intent(in), optional :: vmold + end subroutine mld_d_diag_solver_bld + end interface mld_d_diag_solver_bld + + contains - subroutine d_diag_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) - use psb_base_mod - type(psb_desc_type), intent(in) :: desc_data - class(mld_d_diag_solver_type), intent(inout) :: sv - type(psb_d_vect_type), intent(inout) :: x - type(psb_d_vect_type), intent(inout) :: y - real(psb_dpk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - real(psb_dpk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - - integer :: n_row,n_col - real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act - character :: trans_ - character(len=20) :: name='d_diag_solver_apply' - - call psb_erractionsave(err_act) - - info = psb_success_ - - trans_ = psb_toupper(trans) - select case(trans_) - case('N') - case('T','C') - case default - call psb_errpush(psb_err_iarg_invalid_i_,name) - goto 9999 - end select - - n_row = desc_data%get_local_rows() - n_col = desc_data%get_local_cols() - if (x%get_nrows() < n_row) then - info = 36 - call psb_errpush(info,name,i_err=(/2,nrow,0,0,0/)) - goto 9999 - end if - if (y%get_nrows() < n_row) then - info = 36 - call psb_errpush(info,name,i_err=(/3,nrow,0,0,0/)) - goto 9999 - end if - if (.not.allocated(sv%dv)) then - info = 1124 - call psb_errpush(info,name,a_err="preconditioner: D") - goto 9999 - end if - if (sv%dv%get_nrows() < n_row) then - info = 1124 - call psb_errpush(info,name,a_err="preconditioner: D") - goto 9999 - end if - - call y%mlt(alpha,sv%dv,x,beta,info) - - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='vect%mlt') - goto 9999 - end if - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - - end subroutine d_diag_solver_apply_vect - - subroutine d_diag_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) - use psb_base_mod - type(psb_desc_type), intent(in) :: desc_data - class(mld_d_diag_solver_type), intent(in) :: sv - real(psb_dpk_), intent(inout) :: x(:) - real(psb_dpk_), intent(inout) :: y(:) - real(psb_dpk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - real(psb_dpk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - - integer :: n_row,n_col - real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act - character :: trans_ - character(len=20) :: name='d_diag_solver_apply' - - call psb_erractionsave(err_act) - - info = psb_success_ - - trans_ = psb_toupper(trans) - select case(trans_) - case('N') - case('T','C') - case default - call psb_errpush(psb_err_iarg_invalid_i_,name) - goto 9999 - end select - - n_row = desc_data%get_local_rows() - n_col = desc_data%get_local_cols() - - if (beta == dzero) then - - if (alpha == dzero) then - y(1:n_row) = dzero - else if (alpha == done) then - do i=1, n_row - y(i) = sv%d(i) * x(i) - end do - else if (alpha == -done) then - do i=1, n_row - y(i) = -sv%d(i) * x(i) - end do - else - do i=1, n_row - y(i) = alpha * sv%d(i) * x(i) - end do - end if - - else if (beta == done) then - - if (alpha == dzero) then - !y(1:n_row) = dzero - else if (alpha == done) then - do i=1, n_row - y(i) = sv%d(i) * x(i) + y(i) - end do - else if (alpha == -done) then - do i=1, n_row - y(i) = -sv%d(i) * x(i) + y(i) - end do - else - do i=1, n_row - y(i) = alpha * sv%d(i) * x(i) + y(i) - end do - end if - - else if (beta == -done) then - - if (alpha == dzero) then - y(1:n_row) = -y(1:n_row) - else if (alpha == done) then - do i=1, n_row - y(i) = sv%d(i) * x(i) - y(i) - end do - else if (alpha == -done) then - do i=1, n_row - y(i) = -sv%d(i) * x(i) - y(i) - end do - else - do i=1, n_row - y(i) = alpha * sv%d(i) * x(i) - y(i) - end do - end if - - else - - if (alpha == dzero) then - y(1:n_row) = beta *y(1:n_row) - else if (alpha == done) then - do i=1, n_row - y(i) = sv%d(i) * x(i) + beta*y(i) - end do - else if (alpha == -done) then - do i=1, n_row - y(i) = -sv%d(i) * x(i) + beta*y(i) - end do - else - do i=1, n_row - y(i) = alpha * sv%d(i) * x(i) + beta*y(i) - end do - end if - - end if - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - - end subroutine d_diag_solver_apply - - subroutine d_diag_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold) - - use psb_base_mod - - Implicit None - - ! Arguments - type(psb_dspmat_type), intent(in), target :: a - Type(psb_desc_type), Intent(in) :: desc_a - class(mld_d_diag_solver_type), intent(inout) :: sv - character, intent(in) :: upd - integer, intent(out) :: info - type(psb_dspmat_type), intent(in), target, optional :: b - class(psb_d_base_sparse_mat), intent(in), optional :: amold - class(psb_d_base_vect_type), intent(in), optional :: vmold - ! Local variables - integer :: n_row,n_col, nrow_a, nztota - real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act, debug_unit, debug_level - character(len=20) :: name='d_diag_solver_bld', ch_err - - info=psb_success_ - call psb_erractionsave(err_act) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - ictxt = desc_a%get_context() - call psb_info(ictxt, me, np) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),' start' - - - n_row = desc_a%get_local_rows() - nrow_a = a%get_nrows() - if (allocated(sv%d)) then - if (size(sv%d) < n_row) then - deallocate(sv%d) - endif - endif - if (.not.allocated(sv%d)) then - allocate(sv%d(n_row),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 - end if - - endif - - call a%get_diag(sv%d,info) - if (present(b)) then - if (info == psb_success_) call b%get_diag(sv%d(nrow_a+1:), info) - end if - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='get_diag') - goto 9999 - end if - - do i=1,n_row - if (sv%d(i) == dzero) then - sv%d(i) = done - else - sv%d(i) = done/sv%d(i) - end if - end do - allocate(sv%dv,stat=info) - if (info == psb_success_) then - if (present(vmold)) then - allocate(sv%dv%v,mold=vmold,stat=info) - else - allocate(psb_d_base_vect_type :: sv%dv%v,stat=info) - end if - end if - if (info == psb_success_) then - call sv%dv%bld(sv%d) - else - call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate sv%dv') - goto 9999 - end if - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),' end' - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - end subroutine d_diag_solver_bld - subroutine d_diag_solver_seti(sv,what,val,info) - use psb_base_mod - Implicit None ! Arguments @@ -375,34 +132,14 @@ contains character(len=20) :: name='d_diag_solver_seti' info = psb_success_ -!!$ call psb_erractionsave(err_act) -!!$ -!!$ select case(what) -!!$ case(mld_sub_solve_) -!!$ sv%fact_type = val -!!$ case(mld_sub_fillin_) -!!$ sv%fill_in = val -!!$ case default -!!$ write(0,*) name,': Error: invalid WHAT' -!!$ info = -2 -!!$ end select -!!$ -!!$ call psb_erractionrestore(err_act) - return -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if + return + end subroutine d_diag_solver_seti subroutine d_diag_solver_setc(sv,what,val,info) - use psb_base_mod - Implicit None ! Arguments @@ -414,33 +151,12 @@ contains character(len=20) :: name='d_diag_solver_setc' info = psb_success_ -!!$ call psb_erractionsave(err_act) -!!$ -!!$ -!!$ call mld_stringval(val,ival,info) -!!$ if (info == psb_success_) call sv%set(what,ival,info) -!!$ if (info /= psb_success_) then -!!$ info = psb_err_from_subroutine_ -!!$ call psb_errpush(info, name) -!!$ goto 9999 -!!$ end if -!!$ -!!$ call psb_erractionrestore(err_act) - return -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if return end subroutine d_diag_solver_setc subroutine d_diag_solver_setr(sv,what,val,info) - use psb_base_mod - Implicit None ! Arguments @@ -452,33 +168,13 @@ contains character(len=20) :: name='d_diag_solver_setr' info = psb_success_ -!!$ call psb_erractionsave(err_act) -!!$ -!!$ select case(what) -!!$ case(mld_sub_iluthrs_) -!!$ sv%thresh = val -!!$ case default -!!$ write(0,*) name,': Error: invalid WHAT' -!!$ info = -2 -!!$ goto 9999 -!!$ end select -!!$ -!!$ call psb_erractionrestore(err_act) - return -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if return + end subroutine d_diag_solver_setr subroutine d_diag_solver_free(sv,info) - use psb_base_mod - Implicit None ! Arguments @@ -515,8 +211,6 @@ contains subroutine d_diag_solver_descr(sv,info,iout,coarse) - use psb_base_mod - Implicit None ! Arguments @@ -545,7 +239,6 @@ contains end subroutine d_diag_solver_descr function d_diag_solver_sizeof(sv) result(val) - use psb_base_mod implicit none ! Arguments class(mld_d_diag_solver_type), intent(in) :: sv @@ -559,7 +252,6 @@ contains end function d_diag_solver_sizeof function d_diag_solver_get_nzeros(sv) result(val) - use psb_base_mod implicit none ! Arguments class(mld_d_diag_solver_type), intent(in) :: sv diff --git a/mlprec/mld_d_id_solver.f90 b/mlprec/mld_d_id_solver.f90 index c5ac030d..b6cde41a 100644 --- a/mlprec/mld_d_id_solver.f90 +++ b/mlprec/mld_d_id_solver.f90 @@ -50,8 +50,8 @@ module mld_d_id_solver type, extends(mld_d_base_solver_type) :: mld_d_id_solver_type contains procedure, pass(sv) :: build => d_id_solver_bld - procedure, pass(sv) :: apply_v => d_id_solver_apply_vect - procedure, pass(sv) :: apply_a => d_id_solver_apply + procedure, pass(sv) :: apply_v => mld_d_id_solver_apply_vect + procedure, pass(sv) :: apply_a => mld_d_id_solver_apply procedure, pass(sv) :: free => d_id_solver_free procedure, pass(sv) :: seti => d_id_solver_seti procedure, pass(sv) :: setc => d_id_solver_setc @@ -60,111 +60,46 @@ module mld_d_id_solver end type mld_d_id_solver_type - private :: d_id_solver_bld, d_id_solver_apply, & + private :: d_id_solver_bld, & & d_id_solver_free, d_id_solver_seti, & & d_id_solver_setc, d_id_solver_setr,& & d_id_solver_descr + interface mld_d_id_solver_apply_vect + subroutine mld_d_id_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) + import :: psb_desc_type, psb_dspmat_type, psb_d_base_sparse_mat, & + & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, mld_d_id_solver_type + type(psb_desc_type), intent(in) :: desc_data + class(mld_d_id_solver_type), intent(inout) :: sv + type(psb_d_vect_type),intent(inout) :: x + type(psb_d_vect_type),intent(inout) :: y + real(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + end subroutine mld_d_id_solver_apply_vect + end interface mld_d_id_solver_apply_vect + + interface mld_d_id_solver_apply + subroutine mld_d_id_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) + import :: psb_desc_type, psb_dspmat_type, psb_d_base_sparse_mat, & + & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, mld_d_id_solver_type + type(psb_desc_type), intent(in) :: desc_data + class(mld_d_id_solver_type), intent(in) :: sv + real(psb_dpk_),intent(inout) :: x(:) + real(psb_dpk_),intent(inout) :: y(:) + real(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + end subroutine mld_d_id_solver_apply + end interface mld_d_id_solver_apply contains - subroutine d_id_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) - use psb_base_mod - type(psb_desc_type), intent(in) :: desc_data - class(mld_d_id_solver_type), intent(inout) :: sv - type(psb_d_vect_type),intent(inout) :: x - type(psb_d_vect_type),intent(inout) :: y - real(psb_dpk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - real(psb_dpk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - - integer :: n_row,n_col - real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act - character :: trans_ - character(len=20) :: name='d_id_solver_apply' - - call psb_erractionsave(err_act) - - info = psb_success_ - - trans_ = psb_toupper(trans) - select case(trans_) - case('N') - case('T') - case('C') - case default - call psb_errpush(psb_err_iarg_invalid_i_,name) - goto 9999 - end select - - call psb_geaxpby(alpha,x,beta,y,desc_data,info) - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - - end subroutine d_id_solver_apply_vect - - - subroutine d_id_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) - use psb_base_mod - type(psb_desc_type), intent(in) :: desc_data - class(mld_d_id_solver_type), intent(in) :: sv - real(psb_dpk_),intent(inout) :: x(:) - real(psb_dpk_),intent(inout) :: y(:) - real(psb_dpk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - real(psb_dpk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - - integer :: n_row,n_col - real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act - character :: trans_ - character(len=20) :: name='d_id_solver_apply' - - call psb_erractionsave(err_act) - - info = psb_success_ - - trans_ = psb_toupper(trans) - select case(trans_) - case('N') - case('T') - case('C') - case default - call psb_errpush(psb_err_iarg_invalid_i_,name) - goto 9999 - end select - - call psb_geaxpby(alpha,x,beta,y,desc_data,info) - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - - end subroutine d_id_solver_apply subroutine d_id_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold) - use psb_base_mod - Implicit None ! Arguments @@ -183,35 +118,13 @@ contains character(len=20) :: name='d_id_solver_bld', ch_err info=psb_success_ - call psb_erractionsave(err_act) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - ictxt = desc_a%get_context() - call psb_info(ictxt, me, np) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),' start' - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),' end' - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if return end subroutine d_id_solver_bld subroutine d_id_solver_seti(sv,what,val,info) - use psb_base_mod - Implicit None ! Arguments @@ -230,8 +143,6 @@ contains subroutine d_id_solver_setc(sv,what,val,info) - use psb_base_mod - Implicit None ! Arguments @@ -249,8 +160,6 @@ contains subroutine d_id_solver_setr(sv,what,val,info) - use psb_base_mod - Implicit None ! Arguments @@ -269,8 +178,6 @@ contains subroutine d_id_solver_free(sv,info) - use psb_base_mod - Implicit None ! Arguments @@ -286,8 +193,6 @@ contains subroutine d_id_solver_descr(sv,info,iout,coarse) - use psb_base_mod - Implicit None ! Arguments diff --git a/mlprec/mld_s_diag_solver.f90 b/mlprec/mld_s_diag_solver.f90 index 0b57bfaf..8ed27149 100644 --- a/mlprec/mld_s_diag_solver.f90 +++ b/mlprec/mld_s_diag_solver.f90 @@ -51,9 +51,9 @@ module mld_s_diag_solver type(psb_s_vect_type), allocatable :: dv real(psb_spk_), allocatable :: d(:) contains - procedure, pass(sv) :: build => s_diag_solver_bld - procedure, pass(sv) :: apply_v => s_diag_solver_apply_vect - procedure, pass(sv) :: apply_a => s_diag_solver_apply + procedure, pass(sv) :: build => mld_s_diag_solver_bld + procedure, pass(sv) :: apply_v => mld_s_diag_solver_apply_vect + procedure, pass(sv) :: apply_a => mld_s_diag_solver_apply procedure, pass(sv) :: free => s_diag_solver_free procedure, pass(sv) :: seti => s_diag_solver_seti procedure, pass(sv) :: setc => s_diag_solver_setc @@ -64,306 +64,63 @@ module mld_s_diag_solver end type mld_s_diag_solver_type - private :: s_diag_solver_bld, s_diag_solver_apply, & - & s_diag_solver_free, s_diag_solver_seti, & + private :: s_diag_solver_free, s_diag_solver_seti, & & s_diag_solver_setc, s_diag_solver_setr,& & s_diag_solver_descr, s_diag_solver_sizeof,& & s_diag_solver_get_nzeros + interface mld_s_diag_solver_apply_vect + subroutine mld_s_diag_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) + import :: psb_desc_type, psb_sspmat_type, psb_s_base_sparse_mat, & + & psb_s_vect_type, psb_s_base_vect_type, psb_spk_, mld_s_diag_solver_type + type(psb_desc_type), intent(in) :: desc_data + class(mld_s_diag_solver_type), intent(inout) :: sv + type(psb_s_vect_type), intent(inout) :: x + type(psb_s_vect_type), intent(inout) :: y + real(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_spk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + end subroutine mld_s_diag_solver_apply_vect + end interface mld_s_diag_solver_apply_vect + + interface mld_s_diag_solver_apply + subroutine mld_s_diag_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) + import :: psb_desc_type, psb_sspmat_type, psb_s_base_sparse_mat, & + & psb_s_vect_type, psb_s_base_vect_type, psb_spk_, mld_s_diag_solver_type + type(psb_desc_type), intent(in) :: desc_data + class(mld_s_diag_solver_type), intent(in) :: sv + real(psb_spk_), intent(inout) :: x(:) + real(psb_spk_), intent(inout) :: y(:) + real(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_spk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + end subroutine mld_s_diag_solver_apply + end interface mld_s_diag_solver_apply + + interface mld_s_diag_solver_bld + subroutine mld_s_diag_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold) + import :: psb_desc_type, psb_sspmat_type, psb_s_base_sparse_mat, & + & psb_s_vect_type, psb_s_base_vect_type, psb_spk_, mld_s_diag_solver_type + type(psb_sspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_s_diag_solver_type), intent(inout) :: sv + character, intent(in) :: upd + integer, intent(out) :: info + type(psb_sspmat_type), intent(in), target, optional :: b + class(psb_s_base_sparse_mat), intent(in), optional :: amold + class(psb_s_base_vect_type), intent(in), optional :: vmold + end subroutine mld_s_diag_solver_bld + end interface mld_s_diag_solver_bld + + contains - subroutine s_diag_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) - use psb_base_mod - type(psb_desc_type), intent(in) :: desc_data - class(mld_s_diag_solver_type), intent(inout) :: sv - type(psb_s_vect_type), intent(inout) :: x - type(psb_s_vect_type), intent(inout) :: y - real(psb_spk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - real(psb_spk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - - integer :: n_row,n_col - real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act - character :: trans_ - character(len=20) :: name='s_diag_solver_apply' - - call psb_erractionsave(err_act) - - info = psb_success_ - - trans_ = psb_toupper(trans) - select case(trans_) - case('N') - case('T','C') - case default - call psb_errpush(psb_err_iarg_invalid_i_,name) - goto 9999 - end select - - n_row = desc_data%get_local_rows() - n_col = desc_data%get_local_cols() - if (x%get_nrows() < n_row) then - info = 36 - call psb_errpush(info,name,i_err=(/2,nrow,0,0,0/)) - goto 9999 - end if - if (y%get_nrows() < n_row) then - info = 36 - call psb_errpush(info,name,i_err=(/3,nrow,0,0,0/)) - goto 9999 - end if - if (.not.allocated(sv%dv)) then - info = 1124 - call psb_errpush(info,name,a_err="preconditioner: D") - goto 9999 - end if - if (sv%dv%get_nrows() < n_row) then - info = 1124 - call psb_errpush(info,name,a_err="preconditioner: D") - goto 9999 - end if - - call y%mlt(alpha,sv%dv,x,beta,info) - - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='vect%mlt') - goto 9999 - end if - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - - end subroutine s_diag_solver_apply_vect - - subroutine s_diag_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) - use psb_base_mod - type(psb_desc_type), intent(in) :: desc_data - class(mld_s_diag_solver_type), intent(in) :: sv - real(psb_spk_), intent(inout) :: x(:) - real(psb_spk_), intent(inout) :: y(:) - real(psb_spk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - real(psb_spk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - - integer :: n_row,n_col - real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act - character :: trans_ - character(len=20) :: name='s_diag_solver_apply' - - call psb_erractionsave(err_act) - - info = psb_success_ - - trans_ = psb_toupper(trans) - select case(trans_) - case('N') - case('T','C') - case default - call psb_errpush(psb_err_iarg_invalid_i_,name) - goto 9999 - end select - - n_row = desc_data%get_local_rows() - n_col = desc_data%get_local_cols() - - if (beta == szero) then - - if (alpha == szero) then - y(1:n_row) = szero - else if (alpha == sone) then - do i=1, n_row - y(i) = sv%d(i) * x(i) - end do - else if (alpha == -sone) then - do i=1, n_row - y(i) = -sv%d(i) * x(i) - end do - else - do i=1, n_row - y(i) = alpha * sv%d(i) * x(i) - end do - end if - - else if (beta == sone) then - - if (alpha == szero) then - !y(1:n_row) = szero - else if (alpha == sone) then - do i=1, n_row - y(i) = sv%d(i) * x(i) + y(i) - end do - else if (alpha == -sone) then - do i=1, n_row - y(i) = -sv%d(i) * x(i) + y(i) - end do - else - do i=1, n_row - y(i) = alpha * sv%d(i) * x(i) + y(i) - end do - end if - - else if (beta == -sone) then - - if (alpha == szero) then - y(1:n_row) = -y(1:n_row) - else if (alpha == sone) then - do i=1, n_row - y(i) = sv%d(i) * x(i) - y(i) - end do - else if (alpha == -sone) then - do i=1, n_row - y(i) = -sv%d(i) * x(i) - y(i) - end do - else - do i=1, n_row - y(i) = alpha * sv%d(i) * x(i) - y(i) - end do - end if - - else - - if (alpha == szero) then - y(1:n_row) = beta *y(1:n_row) - else if (alpha == sone) then - do i=1, n_row - y(i) = sv%d(i) * x(i) + beta*y(i) - end do - else if (alpha == -sone) then - do i=1, n_row - y(i) = -sv%d(i) * x(i) + beta*y(i) - end do - else - do i=1, n_row - y(i) = alpha * sv%d(i) * x(i) + beta*y(i) - end do - end if - - end if - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - - end subroutine s_diag_solver_apply - - subroutine s_diag_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold) - - use psb_base_mod - - Implicit None - - ! Arguments - type(psb_sspmat_type), intent(in), target :: a - Type(psb_desc_type), Intent(in) :: desc_a - class(mld_s_diag_solver_type), intent(inout) :: sv - character, intent(in) :: upd - integer, intent(out) :: info - type(psb_sspmat_type), intent(in), target, optional :: b - class(psb_s_base_sparse_mat), intent(in), optional :: amold - class(psb_s_base_vect_type), intent(in), optional :: vmold - ! Local variables - integer :: n_row,n_col, nrow_a, nztota - real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act, debug_unit, debug_level - character(len=20) :: name='s_diag_solver_bld', ch_err - - info=psb_success_ - call psb_erractionsave(err_act) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - ictxt = desc_a%get_context() - call psb_info(ictxt, me, np) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),' start' - - - n_row = desc_a%get_local_rows() - nrow_a = a%get_nrows() - if (allocated(sv%d)) then - if (size(sv%d) < n_row) then - deallocate(sv%d) - endif - endif - if (.not.allocated(sv%d)) then - allocate(sv%d(n_row),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 - end if - - endif - - call a%get_diag(sv%d,info) - if (present(b)) then - if (info == psb_success_) call b%get_diag(sv%d(nrow_a+1:), info) - end if - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='get_diag') - goto 9999 - end if - - do i=1,n_row - if (sv%d(i) == szero) then - sv%d(i) = sone - else - sv%d(i) = sone/sv%d(i) - end if - end do - allocate(sv%dv,stat=info) - if (info == psb_success_) then - if (present(vmold)) then - allocate(sv%dv%v,mold=vmold,stat=info) - else - allocate(psb_s_base_vect_type :: sv%dv%v,stat=info) - end if - end if - if (info == psb_success_) then - call sv%dv%bld(sv%d) - else - call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate sv%dv') - goto 9999 - end if - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),' end' - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - end subroutine s_diag_solver_bld - subroutine s_diag_solver_seti(sv,what,val,info) - use psb_base_mod - Implicit None ! Arguments @@ -375,34 +132,14 @@ contains character(len=20) :: name='s_diag_solver_seti' info = psb_success_ -!!$ call psb_erractionsave(err_act) -!!$ -!!$ select case(what) -!!$ case(mld_sub_solve_) -!!$ sv%fact_type = val -!!$ case(mld_sub_fillin_) -!!$ sv%fill_in = val -!!$ case default -!!$ write(0,*) name,': Error: invalid WHAT' -!!$ info = -2 -!!$ end select -!!$ -!!$ call psb_erractionrestore(err_act) - return -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if + return + end subroutine s_diag_solver_seti subroutine s_diag_solver_setc(sv,what,val,info) - use psb_base_mod - Implicit None ! Arguments @@ -414,33 +151,12 @@ contains character(len=20) :: name='s_diag_solver_setc' info = psb_success_ -!!$ call psb_erractionsave(err_act) -!!$ -!!$ -!!$ call mld_stringval(val,ival,info) -!!$ if (info == psb_success_) call sv%set(what,ival,info) -!!$ if (info /= psb_success_) then -!!$ info = psb_err_from_subroutine_ -!!$ call psb_errpush(info, name) -!!$ goto 9999 -!!$ end if -!!$ -!!$ call psb_erractionrestore(err_act) - return -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if return end subroutine s_diag_solver_setc subroutine s_diag_solver_setr(sv,what,val,info) - use psb_base_mod - Implicit None ! Arguments @@ -452,33 +168,13 @@ contains character(len=20) :: name='s_diag_solver_setr' info = psb_success_ -!!$ call psb_erractionsave(err_act) -!!$ -!!$ select case(what) -!!$ case(mld_sub_iluthrs_) -!!$ sv%thresh = val -!!$ case default -!!$ write(0,*) name,': Error: invalid WHAT' -!!$ info = -2 -!!$ goto 9999 -!!$ end select -!!$ -!!$ call psb_erractionrestore(err_act) - return -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if return + end subroutine s_diag_solver_setr subroutine s_diag_solver_free(sv,info) - use psb_base_mod - Implicit None ! Arguments @@ -515,8 +211,6 @@ contains subroutine s_diag_solver_descr(sv,info,iout,coarse) - use psb_base_mod - Implicit None ! Arguments @@ -545,7 +239,6 @@ contains end subroutine s_diag_solver_descr function s_diag_solver_sizeof(sv) result(val) - use psb_base_mod implicit none ! Arguments class(mld_s_diag_solver_type), intent(in) :: sv @@ -559,7 +252,6 @@ contains end function s_diag_solver_sizeof function s_diag_solver_get_nzeros(sv) result(val) - use psb_base_mod implicit none ! Arguments class(mld_s_diag_solver_type), intent(in) :: sv diff --git a/mlprec/mld_s_id_solver.f90 b/mlprec/mld_s_id_solver.f90 index 1d488a3b..93776b87 100644 --- a/mlprec/mld_s_id_solver.f90 +++ b/mlprec/mld_s_id_solver.f90 @@ -50,8 +50,8 @@ module mld_s_id_solver type, extends(mld_s_base_solver_type) :: mld_s_id_solver_type contains procedure, pass(sv) :: build => s_id_solver_bld - procedure, pass(sv) :: apply_v => s_id_solver_apply_vect - procedure, pass(sv) :: apply_a => s_id_solver_apply + procedure, pass(sv) :: apply_v => mld_s_id_solver_apply_vect + procedure, pass(sv) :: apply_a => mld_s_id_solver_apply procedure, pass(sv) :: free => s_id_solver_free procedure, pass(sv) :: seti => s_id_solver_seti procedure, pass(sv) :: setc => s_id_solver_setc @@ -60,111 +60,46 @@ module mld_s_id_solver end type mld_s_id_solver_type - private :: s_id_solver_bld, s_id_solver_apply, & + private :: s_id_solver_bld, & & s_id_solver_free, s_id_solver_seti, & & s_id_solver_setc, s_id_solver_setr,& & s_id_solver_descr + interface mld_s_id_solver_apply_vect + subroutine mld_s_id_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) + import :: psb_desc_type, psb_sspmat_type, psb_s_base_sparse_mat, & + & psb_s_vect_type, psb_s_base_vect_type, psb_spk_, mld_s_id_solver_type + type(psb_desc_type), intent(in) :: desc_data + class(mld_s_id_solver_type), intent(inout) :: sv + type(psb_s_vect_type),intent(inout) :: x + type(psb_s_vect_type),intent(inout) :: y + real(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_spk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + end subroutine mld_s_id_solver_apply_vect + end interface mld_s_id_solver_apply_vect + + interface mld_s_id_solver_apply + subroutine mld_s_id_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) + import :: psb_desc_type, psb_sspmat_type, psb_s_base_sparse_mat, & + & psb_s_vect_type, psb_s_base_vect_type, psb_spk_, mld_s_id_solver_type + type(psb_desc_type), intent(in) :: desc_data + class(mld_s_id_solver_type), intent(in) :: sv + real(psb_spk_),intent(inout) :: x(:) + real(psb_spk_),intent(inout) :: y(:) + real(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_spk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + end subroutine mld_s_id_solver_apply + end interface mld_s_id_solver_apply contains - subroutine s_id_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) - use psb_base_mod - type(psb_desc_type), intent(in) :: desc_data - class(mld_s_id_solver_type), intent(inout) :: sv - type(psb_s_vect_type),intent(inout) :: x - type(psb_s_vect_type),intent(inout) :: y - real(psb_spk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - real(psb_spk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - - integer :: n_row,n_col - real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act - character :: trans_ - character(len=20) :: name='s_id_solver_apply' - - call psb_erractionsave(err_act) - - info = psb_success_ - - trans_ = psb_toupper(trans) - select case(trans_) - case('N') - case('T') - case('C') - case default - call psb_errpush(psb_err_iarg_invalid_i_,name) - goto 9999 - end select - - call psb_geaxpby(alpha,x,beta,y,desc_data,info) - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - - end subroutine s_id_solver_apply_vect - - - subroutine s_id_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) - use psb_base_mod - type(psb_desc_type), intent(in) :: desc_data - class(mld_s_id_solver_type), intent(in) :: sv - real(psb_spk_),intent(inout) :: x(:) - real(psb_spk_),intent(inout) :: y(:) - real(psb_spk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - real(psb_spk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - - integer :: n_row,n_col - real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act - character :: trans_ - character(len=20) :: name='s_id_solver_apply' - - call psb_erractionsave(err_act) - - info = psb_success_ - - trans_ = psb_toupper(trans) - select case(trans_) - case('N') - case('T') - case('C') - case default - call psb_errpush(psb_err_iarg_invalid_i_,name) - goto 9999 - end select - - call psb_geaxpby(alpha,x,beta,y,desc_data,info) - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - - end subroutine s_id_solver_apply subroutine s_id_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold) - use psb_base_mod - Implicit None ! Arguments @@ -183,35 +118,13 @@ contains character(len=20) :: name='s_id_solver_bld', ch_err info=psb_success_ - call psb_erractionsave(err_act) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - ictxt = desc_a%get_context() - call psb_info(ictxt, me, np) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),' start' - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),' end' - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if return end subroutine s_id_solver_bld subroutine s_id_solver_seti(sv,what,val,info) - use psb_base_mod - Implicit None ! Arguments @@ -230,8 +143,6 @@ contains subroutine s_id_solver_setc(sv,what,val,info) - use psb_base_mod - Implicit None ! Arguments @@ -249,8 +160,6 @@ contains subroutine s_id_solver_setr(sv,what,val,info) - use psb_base_mod - Implicit None ! Arguments @@ -269,8 +178,6 @@ contains subroutine s_id_solver_free(sv,info) - use psb_base_mod - Implicit None ! Arguments @@ -286,8 +193,6 @@ contains subroutine s_id_solver_descr(sv,info,iout,coarse) - use psb_base_mod - Implicit None ! Arguments diff --git a/mlprec/mld_z_diag_solver.f90 b/mlprec/mld_z_diag_solver.f90 index 0671e46b..afd2d350 100644 --- a/mlprec/mld_z_diag_solver.f90 +++ b/mlprec/mld_z_diag_solver.f90 @@ -51,9 +51,9 @@ module mld_z_diag_solver type(psb_z_vect_type), allocatable :: dv complex(psb_dpk_), allocatable :: d(:) contains - procedure, pass(sv) :: build => z_diag_solver_bld - procedure, pass(sv) :: apply_v => z_diag_solver_apply_vect - procedure, pass(sv) :: apply_a => z_diag_solver_apply + procedure, pass(sv) :: build => mld_z_diag_solver_bld + procedure, pass(sv) :: apply_v => mld_z_diag_solver_apply_vect + procedure, pass(sv) :: apply_a => mld_z_diag_solver_apply procedure, pass(sv) :: free => z_diag_solver_free procedure, pass(sv) :: seti => z_diag_solver_seti procedure, pass(sv) :: setc => z_diag_solver_setc @@ -64,306 +64,63 @@ module mld_z_diag_solver end type mld_z_diag_solver_type - private :: z_diag_solver_bld, z_diag_solver_apply, & - & z_diag_solver_free, z_diag_solver_seti, & + private :: z_diag_solver_free, z_diag_solver_seti, & & z_diag_solver_setc, z_diag_solver_setr,& & z_diag_solver_descr, z_diag_solver_sizeof,& & z_diag_solver_get_nzeros + interface mld_z_diag_solver_apply_vect + subroutine mld_z_diag_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) + import :: psb_desc_type, psb_zspmat_type, psb_z_base_sparse_mat, & + & psb_z_vect_type, psb_z_base_vect_type, psb_dpk_, mld_z_diag_solver_type + type(psb_desc_type), intent(in) :: desc_data + class(mld_z_diag_solver_type), intent(inout) :: sv + type(psb_z_vect_type), intent(inout) :: x + type(psb_z_vect_type), intent(inout) :: y + complex(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + end subroutine mld_z_diag_solver_apply_vect + end interface mld_z_diag_solver_apply_vect + + interface mld_z_diag_solver_apply + subroutine mld_z_diag_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) + import :: psb_desc_type, psb_zspmat_type, psb_z_base_sparse_mat, & + & psb_z_vect_type, psb_z_base_vect_type, psb_dpk_, mld_z_diag_solver_type + type(psb_desc_type), intent(in) :: desc_data + class(mld_z_diag_solver_type), intent(in) :: sv + complex(psb_dpk_), intent(inout) :: x(:) + complex(psb_dpk_), intent(inout) :: y(:) + complex(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + end subroutine mld_z_diag_solver_apply + end interface mld_z_diag_solver_apply + + interface mld_z_diag_solver_bld + subroutine mld_z_diag_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold) + import :: psb_desc_type, psb_zspmat_type, psb_z_base_sparse_mat, & + & psb_z_vect_type, psb_z_base_vect_type, psb_dpk_, mld_z_diag_solver_type + type(psb_zspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_z_diag_solver_type), intent(inout) :: sv + character, intent(in) :: upd + integer, intent(out) :: info + type(psb_zspmat_type), intent(in), target, optional :: b + class(psb_z_base_sparse_mat), intent(in), optional :: amold + class(psb_z_base_vect_type), intent(in), optional :: vmold + end subroutine mld_z_diag_solver_bld + end interface mld_z_diag_solver_bld + + contains - subroutine z_diag_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) - use psb_base_mod - type(psb_desc_type), intent(in) :: desc_data - class(mld_z_diag_solver_type), intent(inout) :: sv - type(psb_z_vect_type), intent(inout) :: x - type(psb_z_vect_type), intent(inout) :: y - complex(psb_dpk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - complex(psb_dpk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - - integer :: n_row,n_col - complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act - character :: trans_ - character(len=20) :: name='z_diag_solver_apply' - - call psb_erractionsave(err_act) - - info = psb_success_ - - trans_ = psb_toupper(trans) - select case(trans_) - case('N') - case('T','C') - case default - call psb_errpush(psb_err_iarg_invalid_i_,name) - goto 9999 - end select - - n_row = desc_data%get_local_rows() - n_col = desc_data%get_local_cols() - if (x%get_nrows() < n_row) then - info = 36 - call psb_errpush(info,name,i_err=(/2,nrow,0,0,0/)) - goto 9999 - end if - if (y%get_nrows() < n_row) then - info = 36 - call psb_errpush(info,name,i_err=(/3,nrow,0,0,0/)) - goto 9999 - end if - if (.not.allocated(sv%dv)) then - info = 1124 - call psb_errpush(info,name,a_err="preconditioner: D") - goto 9999 - end if - if (sv%dv%get_nrows() < n_row) then - info = 1124 - call psb_errpush(info,name,a_err="preconditioner: D") - goto 9999 - end if - - call y%mlt(alpha,sv%dv,x,beta,info) - - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='vect%mlt') - goto 9999 - end if - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - - end subroutine z_diag_solver_apply_vect - - subroutine z_diag_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) - use psb_base_mod - type(psb_desc_type), intent(in) :: desc_data - class(mld_z_diag_solver_type), intent(in) :: sv - complex(psb_dpk_), intent(inout) :: x(:) - complex(psb_dpk_), intent(inout) :: y(:) - complex(psb_dpk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - complex(psb_dpk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - - integer :: n_row,n_col - complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act - character :: trans_ - character(len=20) :: name='z_diag_solver_apply' - - call psb_erractionsave(err_act) - - info = psb_success_ - - trans_ = psb_toupper(trans) - select case(trans_) - case('N') - case('T','C') - case default - call psb_errpush(psb_err_iarg_invalid_i_,name) - goto 9999 - end select - - n_row = desc_data%get_local_rows() - n_col = desc_data%get_local_cols() - - if (beta == zzero) then - - if (alpha == zzero) then - y(1:n_row) = zzero - else if (alpha == zone) then - do i=1, n_row - y(i) = sv%d(i) * x(i) - end do - else if (alpha == -zone) then - do i=1, n_row - y(i) = -sv%d(i) * x(i) - end do - else - do i=1, n_row - y(i) = alpha * sv%d(i) * x(i) - end do - end if - - else if (beta == zone) then - - if (alpha == zzero) then - !y(1:n_row) = zzero - else if (alpha == zone) then - do i=1, n_row - y(i) = sv%d(i) * x(i) + y(i) - end do - else if (alpha == -zone) then - do i=1, n_row - y(i) = -sv%d(i) * x(i) + y(i) - end do - else - do i=1, n_row - y(i) = alpha * sv%d(i) * x(i) + y(i) - end do - end if - - else if (beta == -zone) then - - if (alpha == zzero) then - y(1:n_row) = -y(1:n_row) - else if (alpha == zone) then - do i=1, n_row - y(i) = sv%d(i) * x(i) - y(i) - end do - else if (alpha == -zone) then - do i=1, n_row - y(i) = -sv%d(i) * x(i) - y(i) - end do - else - do i=1, n_row - y(i) = alpha * sv%d(i) * x(i) - y(i) - end do - end if - - else - - if (alpha == zzero) then - y(1:n_row) = beta *y(1:n_row) - else if (alpha == zone) then - do i=1, n_row - y(i) = sv%d(i) * x(i) + beta*y(i) - end do - else if (alpha == -zone) then - do i=1, n_row - y(i) = -sv%d(i) * x(i) + beta*y(i) - end do - else - do i=1, n_row - y(i) = alpha * sv%d(i) * x(i) + beta*y(i) - end do - end if - - end if - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - - end subroutine z_diag_solver_apply - - subroutine z_diag_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold) - - use psb_base_mod - - Implicit None - - ! Arguments - type(psb_zspmat_type), intent(in), target :: a - Type(psb_desc_type), Intent(in) :: desc_a - class(mld_z_diag_solver_type), intent(inout) :: sv - character, intent(in) :: upd - integer, intent(out) :: info - type(psb_zspmat_type), intent(in), target, optional :: b - class(psb_z_base_sparse_mat), intent(in), optional :: amold - class(psb_z_base_vect_type), intent(in), optional :: vmold - ! Local variables - integer :: n_row,n_col, nrow_a, nztota - complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act, debug_unit, debug_level - character(len=20) :: name='z_diag_solver_bld', ch_err - - info=psb_success_ - call psb_erractionsave(err_act) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - ictxt = desc_a%get_context() - call psb_info(ictxt, me, np) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),' start' - - - n_row = desc_a%get_local_rows() - nrow_a = a%get_nrows() - if (allocated(sv%d)) then - if (size(sv%d) < n_row) then - deallocate(sv%d) - endif - endif - if (.not.allocated(sv%d)) then - allocate(sv%d(n_row),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 - end if - - endif - - call a%get_diag(sv%d,info) - if (present(b)) then - if (info == psb_success_) call b%get_diag(sv%d(nrow_a+1:), info) - end if - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='get_diag') - goto 9999 - end if - - do i=1,n_row - if (sv%d(i) == zzero) then - sv%d(i) = zone - else - sv%d(i) = zone/sv%d(i) - end if - end do - allocate(sv%dv,stat=info) - if (info == psb_success_) then - if (present(vmold)) then - allocate(sv%dv%v,mold=vmold,stat=info) - else - allocate(psb_z_base_vect_type :: sv%dv%v,stat=info) - end if - end if - if (info == psb_success_) then - call sv%dv%bld(sv%d) - else - call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate sv%dv') - goto 9999 - end if - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),' end' - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - end subroutine z_diag_solver_bld - subroutine z_diag_solver_seti(sv,what,val,info) - use psb_base_mod - Implicit None ! Arguments @@ -375,34 +132,14 @@ contains character(len=20) :: name='z_diag_solver_seti' info = psb_success_ -!!$ call psb_erractionsave(err_act) -!!$ -!!$ select case(what) -!!$ case(mld_sub_solve_) -!!$ sv%fact_type = val -!!$ case(mld_sub_fillin_) -!!$ sv%fill_in = val -!!$ case default -!!$ write(0,*) name,': Error: invalid WHAT' -!!$ info = -2 -!!$ end select -!!$ -!!$ call psb_erractionrestore(err_act) - return -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if + return + end subroutine z_diag_solver_seti subroutine z_diag_solver_setc(sv,what,val,info) - use psb_base_mod - Implicit None ! Arguments @@ -414,33 +151,12 @@ contains character(len=20) :: name='z_diag_solver_setc' info = psb_success_ -!!$ call psb_erractionsave(err_act) -!!$ -!!$ -!!$ call mld_stringval(val,ival,info) -!!$ if (info == psb_success_) call sv%set(what,ival,info) -!!$ if (info /= psb_success_) then -!!$ info = psb_err_from_subroutine_ -!!$ call psb_errpush(info, name) -!!$ goto 9999 -!!$ end if -!!$ -!!$ call psb_erractionrestore(err_act) - return -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if return end subroutine z_diag_solver_setc subroutine z_diag_solver_setr(sv,what,val,info) - use psb_base_mod - Implicit None ! Arguments @@ -452,33 +168,13 @@ contains character(len=20) :: name='z_diag_solver_setr' info = psb_success_ -!!$ call psb_erractionsave(err_act) -!!$ -!!$ select case(what) -!!$ case(mld_sub_iluthrs_) -!!$ sv%thresh = val -!!$ case default -!!$ write(0,*) name,': Error: invalid WHAT' -!!$ info = -2 -!!$ goto 9999 -!!$ end select -!!$ -!!$ call psb_erractionrestore(err_act) - return -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if return + end subroutine z_diag_solver_setr subroutine z_diag_solver_free(sv,info) - use psb_base_mod - Implicit None ! Arguments @@ -515,8 +211,6 @@ contains subroutine z_diag_solver_descr(sv,info,iout,coarse) - use psb_base_mod - Implicit None ! Arguments @@ -545,7 +239,6 @@ contains end subroutine z_diag_solver_descr function z_diag_solver_sizeof(sv) result(val) - use psb_base_mod implicit none ! Arguments class(mld_z_diag_solver_type), intent(in) :: sv @@ -559,7 +252,6 @@ contains end function z_diag_solver_sizeof function z_diag_solver_get_nzeros(sv) result(val) - use psb_base_mod implicit none ! Arguments class(mld_z_diag_solver_type), intent(in) :: sv diff --git a/mlprec/mld_z_id_solver.f90 b/mlprec/mld_z_id_solver.f90 index 680f38b9..f634f891 100644 --- a/mlprec/mld_z_id_solver.f90 +++ b/mlprec/mld_z_id_solver.f90 @@ -50,8 +50,8 @@ module mld_z_id_solver type, extends(mld_z_base_solver_type) :: mld_z_id_solver_type contains procedure, pass(sv) :: build => z_id_solver_bld - procedure, pass(sv) :: apply_v => z_id_solver_apply_vect - procedure, pass(sv) :: apply_a => z_id_solver_apply + procedure, pass(sv) :: apply_v => mld_z_id_solver_apply_vect + procedure, pass(sv) :: apply_a => mld_z_id_solver_apply procedure, pass(sv) :: free => z_id_solver_free procedure, pass(sv) :: seti => z_id_solver_seti procedure, pass(sv) :: setc => z_id_solver_setc @@ -60,111 +60,46 @@ module mld_z_id_solver end type mld_z_id_solver_type - private :: z_id_solver_bld, z_id_solver_apply, & + private :: z_id_solver_bld, & & z_id_solver_free, z_id_solver_seti, & & z_id_solver_setc, z_id_solver_setr,& & z_id_solver_descr + interface mld_z_id_solver_apply_vect + subroutine mld_z_id_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) + import :: psb_desc_type, psb_zspmat_type, psb_z_base_sparse_mat, & + & psb_z_vect_type, psb_z_base_vect_type, psb_dpk_, mld_z_id_solver_type + type(psb_desc_type), intent(in) :: desc_data + class(mld_z_id_solver_type), intent(inout) :: sv + type(psb_z_vect_type),intent(inout) :: x + type(psb_z_vect_type),intent(inout) :: y + complex(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + end subroutine mld_z_id_solver_apply_vect + end interface mld_z_id_solver_apply_vect + + interface mld_z_id_solver_apply + subroutine mld_z_id_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) + import :: psb_desc_type, psb_zspmat_type, psb_z_base_sparse_mat, & + & psb_z_vect_type, psb_z_base_vect_type, psb_dpk_, mld_z_id_solver_type + type(psb_desc_type), intent(in) :: desc_data + class(mld_z_id_solver_type), intent(in) :: sv + complex(psb_dpk_),intent(inout) :: x(:) + complex(psb_dpk_),intent(inout) :: y(:) + complex(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + end subroutine mld_z_id_solver_apply + end interface mld_z_id_solver_apply contains - subroutine z_id_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) - use psb_base_mod - type(psb_desc_type), intent(in) :: desc_data - class(mld_z_id_solver_type), intent(inout) :: sv - type(psb_z_vect_type),intent(inout) :: x - type(psb_z_vect_type),intent(inout) :: y - complex(psb_dpk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - complex(psb_dpk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - - integer :: n_row,n_col - complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act - character :: trans_ - character(len=20) :: name='z_id_solver_apply' - - call psb_erractionsave(err_act) - - info = psb_success_ - - trans_ = psb_toupper(trans) - select case(trans_) - case('N') - case('T') - case('C') - case default - call psb_errpush(psb_err_iarg_invalid_i_,name) - goto 9999 - end select - - call psb_geaxpby(alpha,x,beta,y,desc_data,info) - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - - end subroutine z_id_solver_apply_vect - - - subroutine z_id_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) - use psb_base_mod - type(psb_desc_type), intent(in) :: desc_data - class(mld_z_id_solver_type), intent(in) :: sv - complex(psb_dpk_),intent(inout) :: x(:) - complex(psb_dpk_),intent(inout) :: y(:) - complex(psb_dpk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - complex(psb_dpk_),target, intent(inout) :: work(:) - integer, intent(out) :: info - - integer :: n_row,n_col - complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) - integer :: ictxt,np,me,i, err_act - character :: trans_ - character(len=20) :: name='z_id_solver_apply' - - call psb_erractionsave(err_act) - - info = psb_success_ - - trans_ = psb_toupper(trans) - select case(trans_) - case('N') - case('T') - case('C') - case default - call psb_errpush(psb_err_iarg_invalid_i_,name) - goto 9999 - end select - - call psb_geaxpby(alpha,x,beta,y,desc_data,info) - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - - end subroutine z_id_solver_apply subroutine z_id_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold) - use psb_base_mod - Implicit None ! Arguments @@ -183,35 +118,13 @@ contains character(len=20) :: name='z_id_solver_bld', ch_err info=psb_success_ - call psb_erractionsave(err_act) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - ictxt = desc_a%get_context() - call psb_info(ictxt, me, np) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),' start' - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),' end' - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if return end subroutine z_id_solver_bld subroutine z_id_solver_seti(sv,what,val,info) - use psb_base_mod - Implicit None ! Arguments @@ -230,8 +143,6 @@ contains subroutine z_id_solver_setc(sv,what,val,info) - use psb_base_mod - Implicit None ! Arguments @@ -249,8 +160,6 @@ contains subroutine z_id_solver_setr(sv,what,val,info) - use psb_base_mod - Implicit None ! Arguments @@ -269,8 +178,6 @@ contains subroutine z_id_solver_free(sv,info) - use psb_base_mod - Implicit None ! Arguments @@ -286,8 +193,6 @@ contains subroutine z_id_solver_descr(sv,info,iout,coarse) - use psb_base_mod - Implicit None ! Arguments