diff --git a/mlprec/impl/smoother/Makefile b/mlprec/impl/smoother/Makefile index 6201076e..4b671ad2 100644 --- a/mlprec/impl/smoother/Makefile +++ b/mlprec/impl/smoother/Makefile @@ -22,6 +22,10 @@ mld_c_as_smoother_free.o \ mld_c_as_smoother_setc.o \ mld_c_as_smoother_seti.o \ mld_c_as_smoother_setr.o \ +mld_c_as_smoother_prol_a.o \ +mld_c_as_smoother_prol_v.o \ +mld_c_as_smoother_restr_a.o \ +mld_c_as_smoother_restr_v.o \ mld_c_base_smoother_apply.o \ mld_c_base_smoother_apply_vect.o \ mld_c_base_smoother_bld.o \ @@ -57,6 +61,10 @@ mld_d_as_smoother_free.o \ mld_d_as_smoother_setc.o \ mld_d_as_smoother_seti.o \ mld_d_as_smoother_setr.o \ +mld_d_as_smoother_prol_a.o \ +mld_d_as_smoother_prol_v.o \ +mld_d_as_smoother_restr_a.o \ +mld_d_as_smoother_restr_v.o \ mld_d_base_smoother_apply.o \ mld_d_base_smoother_apply_vect.o \ mld_d_base_smoother_bld.o \ @@ -92,6 +100,10 @@ mld_s_as_smoother_free.o \ mld_s_as_smoother_setc.o \ mld_s_as_smoother_seti.o \ mld_s_as_smoother_setr.o \ +mld_s_as_smoother_prol_a.o \ +mld_s_as_smoother_prol_v.o \ +mld_s_as_smoother_restr_a.o \ +mld_s_as_smoother_restr_v.o \ mld_s_base_smoother_apply.o \ mld_s_base_smoother_apply_vect.o \ mld_s_base_smoother_bld.o \ @@ -127,6 +139,10 @@ mld_z_as_smoother_free.o \ mld_z_as_smoother_setc.o \ mld_z_as_smoother_seti.o \ mld_z_as_smoother_setr.o \ +mld_z_as_smoother_prol_a.o \ +mld_z_as_smoother_prol_v.o \ +mld_z_as_smoother_restr_a.o \ +mld_z_as_smoother_restr_v.o \ mld_z_base_smoother_apply.o \ mld_z_base_smoother_apply_vect.o \ mld_z_base_smoother_bld.o \ diff --git a/mlprec/impl/smoother/mld_c_as_smoother_apply.f90 b/mlprec/impl/smoother/mld_c_as_smoother_apply.f90 index 830f8fef..c41b943e 100644 --- a/mlprec/impl/smoother/mld_c_as_smoother_apply.f90 +++ b/mlprec/impl/smoother/mld_c_as_smoother_apply.f90 @@ -54,15 +54,15 @@ subroutine mld_c_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,& complex(psb_spk_),intent(inout), optional :: initu(:) integer(psb_ipk_) :: n_row,n_col, nrow_d, i - complex(psb_spk_), pointer :: ww(:), aux(:) - complex(psb_spk_), allocatable :: tx(:),ty(:) + complex(psb_spk_), pointer :: aux(:) + complex(psb_spk_), allocatable :: tx(:),ty(:), ww(:) integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5) character :: trans_, init_ character(len=20) :: name='c_as_smoother_apply', ch_err call psb_erractionsave(err_act) - info = psb_success_ + info = psb_success_ ictxt = desc_data%get_context() call psb_info(ictxt,me,np) @@ -89,25 +89,14 @@ subroutine mld_c_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,& end if - n_row = sm%desc_data%get_local_rows() - n_col = sm%desc_data%get_local_cols() + n_row = sm%desc_data%get_local_rows() + n_col = sm%desc_data%get_local_cols() nrow_d = desc_data%get_local_rows() isz = max(n_row,N_COL) - if ((6*isz) <= size(work)) then - ww => work(1:isz) - aux => work(3*isz+1:) - else if ((4*isz) <= size(work)) then + if ((4*isz) <= size(work)) then aux => work(1:) - allocate(ww(isz),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_alloc_request_,name,& - & i_err=(/3*isz,izero,izero,izero,izero/),& - & a_err='complex(psb_spk_)') - goto 9999 - end if - else if ((3*isz) <= size(work)) then - ww => work(1:isz) + else allocate(aux(4*isz),stat=info) if (info /= psb_success_) then call psb_errpush(psb_err_alloc_request_,name,& @@ -115,29 +104,11 @@ subroutine mld_c_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,& & a_err='complex(psb_spk_)') goto 9999 end if - else - allocate(ww(isz), aux(4*isz),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_alloc_request_,name,& - & i_err=(/4*isz,izero,izero,izero,izero/),& - & a_err='complex(psb_spk_)') - goto 9999 - end if - endif - if (sweeps == 0) then - + if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then ! - ! K^0 = I - ! zero sweeps of any smoother is just the identity. - ! - call psb_geaxpby(alpha,x,beta,y,desc_data,info) - - else if ((sm%novr == 0).and.(sweeps == 1).and.(.not.sm%sv%is_iterative())) then - ! - ! Shortcut: in this case it's just the same - ! as Block Jacobi. Moreover, if .not.sv%is_iterative, there's no need to pass init + ! Shortcut: in this case there is nothing else to be done. ! call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) @@ -147,380 +118,114 @@ subroutine mld_c_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,& goto 9999 endif - else - + else if (sweeps >= 0) then + ! + ! + ! Apply multiple sweeps of an AS solver + ! to compute an approximate solution of a linear system. + ! + ! + call psb_geasb(tx,sm%desc_data,info) + call psb_geasb(ty,sm%desc_data,info) + call psb_geasb(ww,sm%desc_data,info) - call psb_geasb(tx,desc_data,info) - call psb_geasb(ty,desc_data,info) + ! + ! Unroll the first iteration and fold it inside SELECT CASE + ! this will save one SPMM when INIT=Z, and will be + ! significant when sweeps=1 (a common case) + ! + call psb_geaxpby(cone,x,czero,tx,desc_data,info) + if (info == 0) call sm%apply_restr(tx,trans_,aux,info) + if (info == 0) call psb_geaxpby(cone,tx,czero,ww,sm%desc_data,info) select case (init_) - case('Z') - tx(:) = czero + case('Z') + call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Z') + case('Y') - call psb_geaxpby(cone,y,czero,tx,desc_data,info) + call psb_geaxpby(cone,y,czero,ty,desc_data,info) + if (info == 0) call sm%apply_restr(ty,trans_,aux,info) + if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,& + & work=aux,trans=trans_) + call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,info,init='Y') + case('U') if (.not.present(initu)) then call psb_errpush(psb_err_internal_error_,name,& & a_err='missing initu to smoother_apply') goto 9999 end if - call psb_geaxpby(cone,initu,czero,tx,desc_data,info) + call psb_geaxpby(cone,initu,czero,ty,desc_data,info) + if (info == 0) call sm%apply_restr(ty,trans_,aux,info) + if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,& + & work=aux,trans=trans_) + call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,info,init='Y') + case default call psb_errpush(psb_err_internal_error_,name,& & a_err='wrong init to smoother_apply') goto 9999 end select - - - if (sweeps == 1) then - - select case(trans_) - case('N') - ! - ! Get the overlap entries of tx (tx == x) - ! - if (sm%restr == psb_halo_) then - call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - else if (sm%restr /= psb_none_) then - call psb_errpush(psb_err_internal_error_,name,& - &a_err='Invalid mld_sub_restr_') - goto 9999 - end if - - - case('T','C') - ! - ! With transpose, we have to do it here - ! - - select case (sm%prol) - - case(psb_none_) - ! - ! Do nothing - - case(psb_sum_) - ! - ! The transpose of sum is halo - ! - call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - - case(psb_avg_) - ! - ! Tricky one: first we have to scale the overlap entries, - ! which we can do by assignind mode=0, i.e. no communication - ! (hence only scaling), then we do the halo - ! - call psb_ovrl(tx,sm%desc_data,info,& - & update=psb_avg_,work=aux,mode=izero) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_prol_') - goto 9999 - end select - - - case default - info=psb_err_iarg_invalid_i_ - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - - call sm%sv%apply(cone,tx,czero,ty,sm%desc_data,trans_,aux,info,init='Y') - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in sub_aply Jacobi Sweeps = 1') - goto 9999 - endif - - select case(trans_) - case('N') - - select case (sm%prol) - - case(psb_none_) - ! - ! Would work anyway, but since it is supposed to do nothing ... - ! call psb_ovrl(ty,sm%desc_data,info,& - ! & update=sm%prol,work=aux) - - - case(psb_sum_,psb_avg_) - ! - ! Update the overlap of ty - ! - call psb_ovrl(ty,sm%desc_data,info,& - & update=sm%prol,work=aux) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_prol_') - goto 9999 - end select - - case('T','C') - ! - ! With transpose, we have to do it here - ! - if (sm%restr == psb_halo_) then - call psb_ovrl(ty,sm%desc_data,info,& - & update=psb_sum_,work=aux) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - else if (sm%restr /= psb_none_) then - call psb_errpush(psb_err_internal_error_,name,& - &a_err='Invalid mld_sub_restr_') - goto 9999 - end if - - case default - info=psb_err_iarg_invalid_i_ - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - - - - else if (sweeps > 1) then - - ! - ! - ! Apply prec%iprcparm(mld_smoother_sweeps_) sweeps of a block-Jacobi solver - ! to compute an approximate solution of a linear system. + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in sub_aply Jacobi Sweeps = 1') + goto 9999 + endif + + do i=1, sweeps-1 ! + ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the + ! block diagonal part and the remaining part of the local matrix + ! and Y(j) is the approximate solution at sweep j. ! - select case (init_) - case('Z') - ty = czero - case('Y') - call psb_geaxpby(cone,y,czero,ty,sm%desc_data,info) - case('U') - if (.not.present(initu)) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='missing initu to smoother_apply') - goto 9999 - end if - call psb_geaxpby(cone,initu,czero,ty,sm%desc_data,info) - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='wrong init to smoother_apply') - goto 9999 - end select + if (info == 0) call psb_geaxpby(cone,tx,czero,ww,sm%desc_data,info) + if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,& + & work=aux,trans=trans_) - do i=1, sweeps - select case(trans_) - case('N') - ! - ! Get the overlap entries of tx (tx == x) - ! - if (sm%restr == psb_halo_) then - call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - else if (sm%restr /= psb_none_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_restr_') - goto 9999 - end if - - - case('T','C') - ! - ! With transpose, we have to do it here - ! - - select case (sm%prol) - - case(psb_none_) - ! - ! Do nothing - - case(psb_sum_) - ! - ! The transpose of sum is halo - ! - call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - - case(psb_avg_) - ! - ! Tricky one: first we have to scale the overlap entries, - ! which we can do by assignind mode=0, i.e. no communication - ! (hence only scaling), then we do the halo - ! - call psb_ovrl(tx,sm%desc_data,info,& - & update=psb_avg_,work=aux,mode=izero) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_prol_') - goto 9999 - end select - - - case default - info=psb_err_iarg_invalid_i_ - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - ! - ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the - ! block diagonal part and the remaining part of the local matrix - ! and Y(j) is the approximate solution at sweep j. - ! - ww(1:n_row) = tx(1:n_row) - call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,& - & work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Y') - - if (info /= psb_success_) exit - - - select case(trans_) - case('N') - - select case (sm%prol) - - case(psb_none_) - ! - ! Would work anyway, but since it is supposed to do nothing ... - ! call psb_ovrl(ty,sm%desc_data,info,& - ! & update=sm%prol,work=aux) - + if (info /= psb_success_) exit - case(psb_sum_,psb_avg_) - ! - ! Update the overlap of ty - ! - call psb_ovrl(ty,sm%desc_data,info,& - & update=sm%prol,work=aux) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_prol_') - goto 9999 - end select - - case('T','C') - ! - ! With transpose, we have to do it here - ! - if (sm%restr == psb_halo_) then - call psb_ovrl(ty,sm%desc_data,info,& - & update=psb_sum_,work=aux) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - else if (sm%restr /= psb_none_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_restr_') - goto 9999 - end if - - case default - info=psb_err_iarg_invalid_i_ - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - end do - - if (info /= psb_success_) then - info=psb_err_internal_error_ - call psb_errpush(info,name,a_err='subsolve with Jacobi sweeps > 1') - goto 9999 - end if - - - else + call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Y') + + if (info /= psb_success_) exit + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + + end do - info = psb_err_iarg_neg_ + if (info /= psb_success_) then + info=psb_err_internal_error_ call psb_errpush(info,name,& - & i_err=(/itwo,sweeps,izero,izero,izero/)) - goto 9999 - - + & a_err='subsolve with Jacobi sweeps > 1') + goto 9999 end if - + ! ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) ! call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + + + else + + info = psb_err_iarg_neg_ + call psb_errpush(info,name,& + & i_err=(/itwo,sweeps,izero,izero,izero/)) + goto 9999 - end if - - - if ((6*isz) <= size(work)) then - else if ((4*isz) <= size(work)) then - deallocate(ww,tx,ty) - else if ((3*isz) <= size(work)) then - deallocate(aux) - else - deallocate(ww,aux,tx,ty) endif + + if (.not.(4*isz <= size(work))) then + deallocate(aux,stat=info) + endif + if (info ==0) deallocate(ww,tx,ty,stat=info) + if (info /= 0) then + info = psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + end if + call psb_erractionrestore(err_act) return diff --git a/mlprec/impl/smoother/mld_c_as_smoother_apply_vect.f90 b/mlprec/impl/smoother/mld_c_as_smoother_apply_vect.f90 index d3fcd8a9..4ea0dede 100644 --- a/mlprec/impl/smoother/mld_c_as_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_c_as_smoother_apply_vect.f90 @@ -54,12 +54,11 @@ subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& type(psb_c_vect_type),intent(inout), optional :: initu integer(psb_ipk_) :: n_row,n_col, nrow_d, i - complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) - complex(psb_spk_), allocatable :: vx(:) - type(psb_c_vect_type) :: vtx, vty, vww + complex(psb_spk_), pointer :: aux(:) + type(psb_c_vect_type) :: tx, ty, ww integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5) character :: trans_, init_ - character(len=20) :: name='c_as_smoother_apply', ch_err + character(len=20) :: name='c_as_smoother_apply_v', ch_err call psb_erractionsave(err_act) @@ -95,55 +94,21 @@ subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& nrow_d = desc_data%get_local_rows() isz = max(n_row,N_COL) - if ((6*isz) <= size(work)) then - ww => work(1:isz) - tx => work(isz+1:2*isz) - ty => work(2*isz+1:3*isz) - aux => work(3*isz+1:) - else if ((4*isz) <= size(work)) then - aux => work(1:) - allocate(ww(isz),tx(isz),ty(isz),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_alloc_request_,name,& - & i_err=(/3*isz,izero,izero,izero,izero/),& - & a_err='complex(psb_spk_)') - goto 9999 - end if - else if ((3*isz) <= size(work)) then - ww => work(1:isz) - tx => work(isz+1:2*isz) - ty => work(2*isz+1:3*isz) - allocate(aux(4*isz),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_alloc_request_,name,& - & i_err=(/4*isz,izero,izero,izero,izero/),& - & a_err='complex(psb_spk_)') - goto 9999 - end if + if (4*isz <= size(work)) then + aux => work(:) else - allocate(ww(isz),tx(isz),ty(isz),& - &aux(4*isz),stat=info) + allocate(aux(4*isz),stat=info) if (info /= psb_success_) then call psb_errpush(psb_err_alloc_request_,name,& & i_err=(/4*isz,izero,izero,izero,izero/),& & a_err='complex(psb_spk_)') goto 9999 end if - endif - if (sweeps == 0) then - + if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then ! - ! K^0 = I - ! zero sweeps of any smoother is just the identity. - ! - call psb_geaxpby(alpha,x,beta,y,desc_data,info) - - else if ((sm%novr == 0).and.(sweeps == 1).and.(.not.sm%sv%is_iterative())) then - ! - ! Shortcut: in this case it's just the same - ! as Block Jacobi. Moreover, if .not.sv%is_iterative, there's no need to pass init + ! Shortcut: in this case there is nothing else to be done. ! call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) @@ -153,388 +118,121 @@ subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& goto 9999 endif - else + else if (sweeps >= 0) then + ! + ! + ! Apply multiple sweeps of an AS solver + ! to compute an approximate solution of a linear system. + ! + ! + call psb_geasb(tx,sm%desc_data,info,mold=x%v,scratch=.true.) + call psb_geasb(ty,sm%desc_data,info,mold=x%v,scratch=.true.) + call psb_geasb(ww,sm%desc_data,info,mold=x%v,scratch=.true.) - call psb_geasb(vtx,sm%desc_data,info,mold=x%v,scratch=.true.) - call psb_geasb(vty,sm%desc_data,info,mold=x%v,scratch=.true.) - call psb_geasb(vww,sm%desc_data,info,mold=x%v,scratch=.true.) + ! + ! Unroll the first iteration and fold it inside SELECT CASE + ! this will save one SPMM when INIT=Z, and will be + ! significant when sweeps=1 (a common case) + ! + call psb_geaxpby(cone,x,czero,tx,desc_data,info) + if (info == 0) call sm%apply_restr(tx,trans_,aux,info) + if (info == 0) call psb_geaxpby(cone,tx,czero,ww,sm%desc_data,info) select case (init_) - case('Z') - call vtx%zero() + case('Z') + call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Z') + case('Y') - call psb_geaxpby(cone,y,czero,vtx,desc_data,info) + call psb_geaxpby(cone,y,czero,ty,desc_data,info) + if (info == 0) call sm%apply_restr(ty,trans_,aux,info) + if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,& + & work=aux,trans=trans_) + call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,info,init='Y') + case('U') if (.not.present(initu)) then call psb_errpush(psb_err_internal_error_,name,& & a_err='missing initu to smoother_apply') goto 9999 end if - call psb_geaxpby(cone,initu,czero,vtx,desc_data,info) + call psb_geaxpby(cone,initu,czero,ty,desc_data,info) + if (info == 0) call sm%apply_restr(ty,trans_,aux,info) + if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,& + & work=aux,trans=trans_) + call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,info,init='Y') + case default call psb_errpush(psb_err_internal_error_,name,& & a_err='wrong init to smoother_apply') goto 9999 end select - - if (sweeps == 1) then - - select case(trans_) - case('N') - ! - ! Get the overlap entries of tx (tx == x) - ! - if (sm%restr == psb_halo_) then - call psb_halo(vtx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - else if (sm%restr /= psb_none_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_restr_') - goto 9999 - end if - - - case('T','C') - ! - ! With transpose, we have to do it here - ! - - select case (sm%prol) - - case(psb_none_) - ! - ! Do nothing - - case(psb_sum_) - ! - ! The transpose of sum is halo - ! - call psb_halo(vtx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - - case(psb_avg_) - ! - ! Tricky one: first we have to scale the overlap entries, - ! which we can do by assignind mode=0, i.e. no communication - ! (hence only scaling), then we do the halo - ! - call psb_ovrl(vtx,sm%desc_data,info,& - & update=psb_avg_,work=aux,mode=izero) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - call psb_halo(vtx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_prol_') - goto 9999 - end select - - - case default - info=psb_err_iarg_invalid_i_ - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - - call sm%sv%apply(cone,vtx,czero,vty,sm%desc_data,trans_,aux,info,init='Y') - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in sub_aply Jacobi Sweeps = 1') - goto 9999 - endif - - select case(trans_) - case('N') - - select case (sm%prol) - - case(psb_none_) - ! - ! Would work anyway, but since it is supposed to do nothing ... - ! call psb_ovrl(ty,sm%desc_data,info,& - ! & update=sm%prol,work=aux) - - - case(psb_sum_,psb_avg_) - ! - ! Update the overlap of ty - ! - call psb_ovrl(vty,sm%desc_data,info,& - & update=sm%prol,work=aux) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_prol_') - goto 9999 - end select - - case('T','C') - ! - ! With transpose, we have to do it here - ! - if (sm%restr == psb_halo_) then - call psb_ovrl(vty,sm%desc_data,info,& - & update=psb_sum_,work=aux) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - else if (sm%restr /= psb_none_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_restr_') - goto 9999 - end if - - case default - info=psb_err_iarg_invalid_i_ - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - - - - else if (sweeps > 1) then - - ! - ! - ! Apply prec%iprcparm(mld_smoother_sweeps_) sweeps of a block-Jacobi solver - ! to compute an approximate solution of a linear system. + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in sub_aply Jacobi Sweeps = 1') + goto 9999 + endif + + do i=1, sweeps-1 ! + ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the + ! block diagonal part and the remaining part of the local matrix + ! and Y(j) is the approximate solution at sweep j. ! - select case (init_) - case('Z') - call vty%zero() - case('Y') - call psb_geaxpby(cone,y,czero,vty,sm%desc_data,info) - case('U') - if (.not.present(initu)) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='missing initu to smoother_apply') - goto 9999 - end if - call psb_geaxpby(cone,initu,czero,vty,sm%desc_data,info) - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='wrong init to smoother_apply') - goto 9999 - end select + if (info == 0) call psb_geaxpby(cone,tx,czero,ww,sm%desc_data,info) + if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,& + & work=aux,trans=trans_) - do i=1, sweeps - select case(trans_) - case('N') - ! - ! Get the overlap entries of tx (tx == x) - ! - if (sm%restr == psb_halo_) then - call psb_halo(vtx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - else if (sm%restr /= psb_none_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_restr_') - goto 9999 - end if - - - case('T','C') - ! - ! With transpose, we have to do it here - ! - - select case (sm%prol) - - case(psb_none_) - ! - ! Do nothing - - case(psb_sum_) - ! - ! The transpose of sum is halo - ! - call psb_halo(vtx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - - case(psb_avg_) - ! - ! Tricky one: first we have to scale the overlap entries, - ! which we can do by assignind mode=0, i.e. no communication - ! (hence only scaling), then we do the halo - ! - call psb_ovrl(vtx,sm%desc_data,info,& - & update=psb_avg_,work=aux,mode=izero) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - call psb_halo(vtx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_prol_') - goto 9999 - end select - - - case default - info=psb_err_iarg_invalid_i_ - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - ! - ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the - ! block diagonal part and the remaining part of the local matrix - ! and Y(j) is the approximate solution at sweep j. - ! - call psb_geaxpby(cone,vtx,czero,vww,sm%desc_data,info) - call psb_spmm(-cone,sm%nd,vty,cone,vww,sm%desc_data,info,& - & work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call sm%sv%apply(cone,vww,czero,vty,sm%desc_data,trans_,aux,info,init='Y') - - if (info /= psb_success_) exit - - - select case(trans_) - case('N') - - select case (sm%prol) - - case(psb_none_) - ! - ! Would work anyway, but since it is supposed to do nothing ... - ! call psb_ovrl(ty,sm%desc_data,info,& - ! & update=sm%prol,work=aux) - - - case(psb_sum_,psb_avg_) - ! - ! Update the overlap of ty - ! - call psb_ovrl(vty,sm%desc_data,info,& - & update=sm%prol,work=aux) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_prol_') - goto 9999 - end select + if (info /= psb_success_) exit - case('T','C') - ! - ! With transpose, we have to do it here - ! - if (sm%restr == psb_halo_) then - call psb_ovrl(vty,sm%desc_data,info,& - & update=psb_sum_,work=aux) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - else if (sm%restr /= psb_none_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_restr_') - goto 9999 - end if - - case default - info=psb_err_iarg_invalid_i_ - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - end do - - if (info /= psb_success_) then - info=psb_err_internal_error_ - call psb_errpush(info,name,& - & a_err='subsolve with Jacobi sweeps > 1') - goto 9999 - end if - - - else + call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Y') + + if (info /= psb_success_) exit + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + + end do - info = psb_err_iarg_neg_ + if (info /= psb_success_) then + info=psb_err_internal_error_ call psb_errpush(info,name,& - & i_err=(/itwo,sweeps,izero,izero,izero/)) - goto 9999 - - + & a_err='subsolve with Jacobi sweeps > 1') + goto 9999 end if - + ! ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) ! - call psb_geaxpby(alpha,vty,beta,y,desc_data,info) - - end if + call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + + + else + + info = psb_err_iarg_neg_ + call psb_errpush(info,name,& + & i_err=(/itwo,sweeps,izero,izero,izero/)) + goto 9999 + endif - if ((6*isz) <= size(work)) then - else if ((4*isz) <= size(work)) then - deallocate(ww,tx,ty) - else if ((3*isz) <= size(work)) then - deallocate(aux) - else - deallocate(ww,aux,tx,ty) + + if (.not.(4*isz <= size(work))) then + deallocate(aux,stat=info) endif - call vww%free(info) - call vtx%free(info) - call vty%free(info) + if (info ==0) call ww%free(info) + if (info ==0) call tx%free(info) + if (info ==0) call ty%free(info) + if (info /= 0) then + info = psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + end if call psb_erractionrestore(err_act) return - + 9999 call psb_error_handler(err_act) - + return end subroutine mld_c_as_smoother_apply_vect diff --git a/mlprec/impl/smoother/mld_c_as_smoother_prol_a.f90 b/mlprec/impl/smoother/mld_c_as_smoother_prol_a.f90 new file mode 100644 index 00000000..fd3e89e4 --- /dev/null +++ b/mlprec/impl/smoother/mld_c_as_smoother_prol_a.f90 @@ -0,0 +1,150 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ 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_as_smoother_prol_a(sm,x,trans,work,info,data) + use psb_base_mod + use mld_c_as_smoother, mld_protect_nam => mld_c_as_smoother_prol_a + implicit none + class(mld_c_as_smoother_type), intent(inout) :: sm + complex(psb_spk_), intent(inout) :: x(:) + character(len=1),intent(in) :: trans + complex(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: data + !Local + integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5), data_ + character :: trans_ + character(len=20) :: name='c_as_smther_prol_a', ch_err + + call psb_erractionsave(err_act) + + info = psb_success_ + ictxt = sm%desc_data%get_context() + call psb_info(ictxt,me,np) + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T') + case('C') + case default + info = psb_err_iarg_invalid_i_ + call psb_errpush(info,name) + goto 9999 + end select + + if (present(data)) then + data_ = data + else + data_ = psb_comm_ext_ + end if + + if (.not.allocated(sm%sv)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + + + select case(trans_) + case('N') + + select case (sm%prol) + + case(psb_none_) + ! + ! Would work anyway, but since it is supposed to do nothing ... + ! call psb_ovrl(x,sm%desc_data,info,& + ! & update=sm%prol,work=work) + + + case(psb_sum_,psb_avg_) + ! + ! Update the overlap of x + ! + call psb_ovrl(x,sm%desc_data,info,& + & update=sm%prol,work=work) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_ovrl' + goto 9999 + end if + + case default + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Invalid mld_sub_prol_') + goto 9999 + end select + + case('T','C') + ! + ! With transpose, we have to do it here + ! + if (sm%restr == psb_halo_) then + call psb_ovrl(x,sm%desc_data,info,& + & update=psb_sum_,work=work) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_ovrl' + goto 9999 + end if + else if (sm%restr /= psb_none_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Invalid mld_sub_restr_') + goto 9999 + end if + + case default + info=psb_err_iarg_invalid_i_ + int_err(1)=6 + ch_err(2:2)=trans + goto 9999 + end select + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_c_as_smoother_prol_a + + diff --git a/mlprec/impl/smoother/mld_c_as_smoother_prol_v.f90 b/mlprec/impl/smoother/mld_c_as_smoother_prol_v.f90 new file mode 100644 index 00000000..5cb0ff0a --- /dev/null +++ b/mlprec/impl/smoother/mld_c_as_smoother_prol_v.f90 @@ -0,0 +1,150 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ 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_as_smoother_prol_v(sm,x,trans,work,info,data) + use psb_base_mod + use mld_c_as_smoother, mld_protect_nam => mld_c_as_smoother_prol_v + implicit none + class(mld_c_as_smoother_type), intent(inout) :: sm + type(psb_c_vect_type),intent(inout) :: x + character(len=1),intent(in) :: trans + complex(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: data + !Local + integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5), data_ + character :: trans_ + character(len=20) :: name='c_as_smther_prol_v', ch_err + + call psb_erractionsave(err_act) + + info = psb_success_ + ictxt = sm%desc_data%get_context() + call psb_info(ictxt,me,np) + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T') + case('C') + case default + info = psb_err_iarg_invalid_i_ + call psb_errpush(info,name) + goto 9999 + end select + + if (present(data)) then + data_ = data + else + data_ = psb_comm_ext_ + end if + + if (.not.allocated(sm%sv)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + + + select case(trans_) + case('N') + + select case (sm%prol) + + case(psb_none_) + ! + ! Would work anyway, but since it is supposed to do nothing ... + ! call psb_ovrl(x,sm%desc_data,info,& + ! & update=sm%prol,work=work) + + + case(psb_sum_,psb_avg_) + ! + ! Update the overlap of x + ! + call psb_ovrl(x,sm%desc_data,info,& + & update=sm%prol,work=work) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_ovrl' + goto 9999 + end if + + case default + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Invalid mld_sub_prol_') + goto 9999 + end select + + case('T','C') + ! + ! With transpose, we have to do it here + ! + if (sm%restr == psb_halo_) then + call psb_ovrl(x,sm%desc_data,info,& + & update=psb_sum_,work=work) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_ovrl' + goto 9999 + end if + else if (sm%restr /= psb_none_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Invalid mld_sub_restr_') + goto 9999 + end if + + case default + info=psb_err_iarg_invalid_i_ + int_err(1)=6 + ch_err(2:2)=trans + goto 9999 + end select + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_c_as_smoother_prol_v + + diff --git a/mlprec/impl/smoother/mld_c_as_smoother_restr_a.f90 b/mlprec/impl/smoother/mld_c_as_smoother_restr_a.f90 new file mode 100644 index 00000000..6d34175c --- /dev/null +++ b/mlprec/impl/smoother/mld_c_as_smoother_restr_a.f90 @@ -0,0 +1,169 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ 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_as_smoother_restr_a(sm,x,trans,work,info,data) + use psb_base_mod + use mld_c_as_smoother, mld_protect_nam => mld_c_as_smoother_restr_a + implicit none + class(mld_c_as_smoother_type), intent(inout) :: sm + complex(psb_spk_), intent(inout) :: x(:) + character(len=1),intent(in) :: trans + complex(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: data + !Local + integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5), data_ + character :: trans_ + character(len=20) :: name='c_as_smther_restr_a', ch_err + + call psb_erractionsave(err_act) + + info = psb_success_ + ictxt = sm%desc_data%get_context() + call psb_info(ictxt,me,np) + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T') + case('C') + case default + info = psb_err_iarg_invalid_i_ + call psb_errpush(info,name) + goto 9999 + end select + + if (present(data)) then + data_ = data + else + data_ = psb_comm_ext_ + end if + + if (.not.allocated(sm%sv)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + + select case(trans_) + case('N') + ! + ! Get the overlap entries x + ! + if (sm%restr == psb_halo_) then + call psb_halo(x,sm%desc_data,info,work=work,data=data_) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_halo' + goto 9999 + end if + else if (sm%restr /= psb_none_) then + call psb_errpush(psb_err_internal_error_,name,& + &a_err='Invalid mld_sub_restr_') + goto 9999 + end if + + + case('T','C') + ! + ! With transpose, we have to do it here + ! + + select case (sm%prol) + + case(psb_none_) + ! + ! Do nothing + + case(psb_sum_) + ! + ! The transpose of sum is halo + ! + call psb_halo(x,sm%desc_data,info,work=work,data=data_) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_halo' + goto 9999 + end if + + case(psb_avg_) + ! + ! Tricky one: first we have to scale the overlap entries, + ! which we can do by assignind mode=0, i.e. no communication + ! (hence only scaling), then we do the halo + ! + call psb_ovrl(x,sm%desc_data,info,& + & update=psb_avg_,work=work,mode=izero) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_ovrl' + goto 9999 + end if + call psb_halo(x,sm%desc_data,info,work=work,data=data_) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_halo' + goto 9999 + end if + + case default + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Invalid mld_sub_prol_') + goto 9999 + end select + + + case default + info=psb_err_iarg_invalid_i_ + int_err(1)=6 + ch_err(2:2)=trans + goto 9999 + end select + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_c_as_smoother_restr_a + + diff --git a/mlprec/impl/smoother/mld_c_as_smoother_restr_v.f90 b/mlprec/impl/smoother/mld_c_as_smoother_restr_v.f90 new file mode 100644 index 00000000..a3588a88 --- /dev/null +++ b/mlprec/impl/smoother/mld_c_as_smoother_restr_v.f90 @@ -0,0 +1,169 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ 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_as_smoother_restr_v(sm,x,trans,work,info,data) + use psb_base_mod + use mld_c_as_smoother, mld_protect_nam => mld_c_as_smoother_restr_v + implicit none + class(mld_c_as_smoother_type), intent(inout) :: sm + type(psb_c_vect_type),intent(inout) :: x + character(len=1),intent(in) :: trans + complex(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: data + !Local + integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5), data_ + character :: trans_ + character(len=20) :: name='c_as_smther_restr_v', ch_err + + call psb_erractionsave(err_act) + + info = psb_success_ + ictxt = sm%desc_data%get_context() + call psb_info(ictxt,me,np) + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T') + case('C') + case default + info = psb_err_iarg_invalid_i_ + call psb_errpush(info,name) + goto 9999 + end select + + if (present(data)) then + data_ = data + else + data_ = psb_comm_ext_ + end if + + if (.not.allocated(sm%sv)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + + select case(trans_) + case('N') + ! + ! Get the overlap entries x + ! + if (sm%restr == psb_halo_) then + call psb_halo(x,sm%desc_data,info,work=work,data=data_) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_halo' + goto 9999 + end if + else if (sm%restr /= psb_none_) then + call psb_errpush(psb_err_internal_error_,name,& + &a_err='Invalid mld_sub_restr_') + goto 9999 + end if + + + case('T','C') + ! + ! With transpose, we have to do it here + ! + + select case (sm%prol) + + case(psb_none_) + ! + ! Do nothing + + case(psb_sum_) + ! + ! The transpose of sum is halo + ! + call psb_halo(x,sm%desc_data,info,work=work,data=data_) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_halo' + goto 9999 + end if + + case(psb_avg_) + ! + ! Tricky one: first we have to scale the overlap entries, + ! which we can do by assignind mode=0, i.e. no communication + ! (hence only scaling), then we do the halo + ! + call psb_ovrl(x,sm%desc_data,info,& + & update=psb_avg_,work=work,mode=izero) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_ovrl' + goto 9999 + end if + call psb_halo(x,sm%desc_data,info,work=work,data=data_) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_halo' + goto 9999 + end if + + case default + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Invalid mld_sub_prol_') + goto 9999 + end select + + + case default + info=psb_err_iarg_invalid_i_ + int_err(1)=6 + ch_err(2:2)=trans + goto 9999 + end select + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_c_as_smoother_restr_v + + diff --git a/mlprec/impl/smoother/mld_d_as_smoother_apply.f90 b/mlprec/impl/smoother/mld_d_as_smoother_apply.f90 index 8acfe8a4..5b08489a 100644 --- a/mlprec/impl/smoother/mld_d_as_smoother_apply.f90 +++ b/mlprec/impl/smoother/mld_d_as_smoother_apply.f90 @@ -54,15 +54,15 @@ subroutine mld_d_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,& real(psb_dpk_),intent(inout), optional :: initu(:) integer(psb_ipk_) :: n_row,n_col, nrow_d, i - real(psb_dpk_), pointer :: ww(:), aux(:) - real(psb_dpk_), allocatable :: tx(:),ty(:) + real(psb_dpk_), pointer :: aux(:) + real(psb_dpk_), allocatable :: tx(:),ty(:), ww(:) integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5) character :: trans_, init_ character(len=20) :: name='d_as_smoother_apply', ch_err call psb_erractionsave(err_act) - info = psb_success_ + info = psb_success_ ictxt = desc_data%get_context() call psb_info(ictxt,me,np) @@ -89,25 +89,14 @@ subroutine mld_d_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,& end if - n_row = sm%desc_data%get_local_rows() - n_col = sm%desc_data%get_local_cols() + n_row = sm%desc_data%get_local_rows() + n_col = sm%desc_data%get_local_cols() nrow_d = desc_data%get_local_rows() isz = max(n_row,N_COL) - if ((6*isz) <= size(work)) then - ww => work(1:isz) - aux => work(3*isz+1:) - else if ((4*isz) <= size(work)) then + if ((4*isz) <= size(work)) then aux => work(1:) - allocate(ww(isz),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_alloc_request_,name,& - & i_err=(/3*isz,izero,izero,izero,izero/),& - & a_err='real(psb_dpk_)') - goto 9999 - end if - else if ((3*isz) <= size(work)) then - ww => work(1:isz) + else allocate(aux(4*isz),stat=info) if (info /= psb_success_) then call psb_errpush(psb_err_alloc_request_,name,& @@ -115,29 +104,11 @@ subroutine mld_d_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,& & a_err='real(psb_dpk_)') goto 9999 end if - else - allocate(ww(isz), aux(4*isz),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_alloc_request_,name,& - & i_err=(/4*isz,izero,izero,izero,izero/),& - & a_err='real(psb_dpk_)') - goto 9999 - end if - endif - if (sweeps == 0) then - + if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then ! - ! K^0 = I - ! zero sweeps of any smoother is just the identity. - ! - call psb_geaxpby(alpha,x,beta,y,desc_data,info) - - else if ((sm%novr == 0).and.(sweeps == 1).and.(.not.sm%sv%is_iterative())) then - ! - ! Shortcut: in this case it's just the same - ! as Block Jacobi. Moreover, if .not.sv%is_iterative, there's no need to pass init + ! Shortcut: in this case there is nothing else to be done. ! call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) @@ -147,380 +118,114 @@ subroutine mld_d_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,& goto 9999 endif - else - + else if (sweeps >= 0) then + ! + ! + ! Apply multiple sweeps of an AS solver + ! to compute an approximate solution of a linear system. + ! + ! + call psb_geasb(tx,sm%desc_data,info) + call psb_geasb(ty,sm%desc_data,info) + call psb_geasb(ww,sm%desc_data,info) - call psb_geasb(tx,desc_data,info) - call psb_geasb(ty,desc_data,info) + ! + ! Unroll the first iteration and fold it inside SELECT CASE + ! this will save one SPMM when INIT=Z, and will be + ! significant when sweeps=1 (a common case) + ! + call psb_geaxpby(done,x,dzero,tx,desc_data,info) + if (info == 0) call sm%apply_restr(tx,trans_,aux,info) + if (info == 0) call psb_geaxpby(done,tx,dzero,ww,sm%desc_data,info) select case (init_) - case('Z') - tx(:) = dzero + case('Z') + call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Z') + case('Y') - call psb_geaxpby(done,y,dzero,tx,desc_data,info) + call psb_geaxpby(done,y,dzero,ty,desc_data,info) + if (info == 0) call sm%apply_restr(ty,trans_,aux,info) + if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,& + & work=aux,trans=trans_) + call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,info,init='Y') + case('U') if (.not.present(initu)) then call psb_errpush(psb_err_internal_error_,name,& & a_err='missing initu to smoother_apply') goto 9999 end if - call psb_geaxpby(done,initu,dzero,tx,desc_data,info) + call psb_geaxpby(done,initu,dzero,ty,desc_data,info) + if (info == 0) call sm%apply_restr(ty,trans_,aux,info) + if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,& + & work=aux,trans=trans_) + call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,info,init='Y') + case default call psb_errpush(psb_err_internal_error_,name,& & a_err='wrong init to smoother_apply') goto 9999 end select - - - if (sweeps == 1) then - - select case(trans_) - case('N') - ! - ! Get the overlap entries of tx (tx == x) - ! - if (sm%restr == psb_halo_) then - call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - else if (sm%restr /= psb_none_) then - call psb_errpush(psb_err_internal_error_,name,& - &a_err='Invalid mld_sub_restr_') - goto 9999 - end if - - - case('T','C') - ! - ! With transpose, we have to do it here - ! - - select case (sm%prol) - - case(psb_none_) - ! - ! Do nothing - - case(psb_sum_) - ! - ! The transpose of sum is halo - ! - call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - - case(psb_avg_) - ! - ! Tricky one: first we have to scale the overlap entries, - ! which we can do by assignind mode=0, i.e. no communication - ! (hence only scaling), then we do the halo - ! - call psb_ovrl(tx,sm%desc_data,info,& - & update=psb_avg_,work=aux,mode=izero) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_prol_') - goto 9999 - end select - - - case default - info=psb_err_iarg_invalid_i_ - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - - call sm%sv%apply(done,tx,dzero,ty,sm%desc_data,trans_,aux,info,init='Y') - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in sub_aply Jacobi Sweeps = 1') - goto 9999 - endif - - select case(trans_) - case('N') - - select case (sm%prol) - - case(psb_none_) - ! - ! Would work anyway, but since it is supposed to do nothing ... - ! call psb_ovrl(ty,sm%desc_data,info,& - ! & update=sm%prol,work=aux) - - - case(psb_sum_,psb_avg_) - ! - ! Update the overlap of ty - ! - call psb_ovrl(ty,sm%desc_data,info,& - & update=sm%prol,work=aux) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_prol_') - goto 9999 - end select - - case('T','C') - ! - ! With transpose, we have to do it here - ! - if (sm%restr == psb_halo_) then - call psb_ovrl(ty,sm%desc_data,info,& - & update=psb_sum_,work=aux) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - else if (sm%restr /= psb_none_) then - call psb_errpush(psb_err_internal_error_,name,& - &a_err='Invalid mld_sub_restr_') - goto 9999 - end if - - case default - info=psb_err_iarg_invalid_i_ - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - - - - else if (sweeps > 1) then - - ! - ! - ! Apply prec%iprcparm(mld_smoother_sweeps_) sweeps of a block-Jacobi solver - ! to compute an approximate solution of a linear system. + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in sub_aply Jacobi Sweeps = 1') + goto 9999 + endif + + do i=1, sweeps-1 ! + ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the + ! block diagonal part and the remaining part of the local matrix + ! and Y(j) is the approximate solution at sweep j. ! - select case (init_) - case('Z') - ty = dzero - case('Y') - call psb_geaxpby(done,y,dzero,ty,sm%desc_data,info) - case('U') - if (.not.present(initu)) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='missing initu to smoother_apply') - goto 9999 - end if - call psb_geaxpby(done,initu,dzero,ty,sm%desc_data,info) - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='wrong init to smoother_apply') - goto 9999 - end select + if (info == 0) call psb_geaxpby(done,tx,dzero,ww,sm%desc_data,info) + if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,& + & work=aux,trans=trans_) - do i=1, sweeps - select case(trans_) - case('N') - ! - ! Get the overlap entries of tx (tx == x) - ! - if (sm%restr == psb_halo_) then - call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - else if (sm%restr /= psb_none_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_restr_') - goto 9999 - end if - - - case('T','C') - ! - ! With transpose, we have to do it here - ! - - select case (sm%prol) - - case(psb_none_) - ! - ! Do nothing - - case(psb_sum_) - ! - ! The transpose of sum is halo - ! - call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - - case(psb_avg_) - ! - ! Tricky one: first we have to scale the overlap entries, - ! which we can do by assignind mode=0, i.e. no communication - ! (hence only scaling), then we do the halo - ! - call psb_ovrl(tx,sm%desc_data,info,& - & update=psb_avg_,work=aux,mode=izero) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_prol_') - goto 9999 - end select - - - case default - info=psb_err_iarg_invalid_i_ - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - ! - ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the - ! block diagonal part and the remaining part of the local matrix - ! and Y(j) is the approximate solution at sweep j. - ! - ww(1:n_row) = tx(1:n_row) - call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,& - & work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Y') - - if (info /= psb_success_) exit - - - select case(trans_) - case('N') - - select case (sm%prol) - - case(psb_none_) - ! - ! Would work anyway, but since it is supposed to do nothing ... - ! call psb_ovrl(ty,sm%desc_data,info,& - ! & update=sm%prol,work=aux) - + if (info /= psb_success_) exit - case(psb_sum_,psb_avg_) - ! - ! Update the overlap of ty - ! - call psb_ovrl(ty,sm%desc_data,info,& - & update=sm%prol,work=aux) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_prol_') - goto 9999 - end select - - case('T','C') - ! - ! With transpose, we have to do it here - ! - if (sm%restr == psb_halo_) then - call psb_ovrl(ty,sm%desc_data,info,& - & update=psb_sum_,work=aux) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - else if (sm%restr /= psb_none_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_restr_') - goto 9999 - end if - - case default - info=psb_err_iarg_invalid_i_ - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - end do - - if (info /= psb_success_) then - info=psb_err_internal_error_ - call psb_errpush(info,name,a_err='subsolve with Jacobi sweeps > 1') - goto 9999 - end if - - - else + call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Y') + + if (info /= psb_success_) exit + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + + end do - info = psb_err_iarg_neg_ + if (info /= psb_success_) then + info=psb_err_internal_error_ call psb_errpush(info,name,& - & i_err=(/itwo,sweeps,izero,izero,izero/)) - goto 9999 - - + & a_err='subsolve with Jacobi sweeps > 1') + goto 9999 end if - + ! ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) ! call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + + + else + + info = psb_err_iarg_neg_ + call psb_errpush(info,name,& + & i_err=(/itwo,sweeps,izero,izero,izero/)) + goto 9999 - end if - - - if ((6*isz) <= size(work)) then - else if ((4*isz) <= size(work)) then - deallocate(ww,tx,ty) - else if ((3*isz) <= size(work)) then - deallocate(aux) - else - deallocate(ww,aux,tx,ty) endif + + if (.not.(4*isz <= size(work))) then + deallocate(aux,stat=info) + endif + if (info ==0) deallocate(ww,tx,ty,stat=info) + if (info /= 0) then + info = psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + end if + call psb_erractionrestore(err_act) return diff --git a/mlprec/impl/smoother/mld_d_as_smoother_apply_vect.f90 b/mlprec/impl/smoother/mld_d_as_smoother_apply_vect.f90 index 58442077..7e6d4473 100644 --- a/mlprec/impl/smoother/mld_d_as_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_d_as_smoother_apply_vect.f90 @@ -54,12 +54,11 @@ subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& type(psb_d_vect_type),intent(inout), optional :: initu integer(psb_ipk_) :: n_row,n_col, nrow_d, i - real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) - real(psb_dpk_), allocatable :: vx(:) - type(psb_d_vect_type) :: vtx, vty, vww + real(psb_dpk_), pointer :: aux(:) + type(psb_d_vect_type) :: tx, ty, ww integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5) character :: trans_, init_ - character(len=20) :: name='d_as_smoother_apply', ch_err + character(len=20) :: name='d_as_smoother_apply_v', ch_err call psb_erractionsave(err_act) @@ -95,55 +94,21 @@ subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& nrow_d = desc_data%get_local_rows() isz = max(n_row,N_COL) - if ((6*isz) <= size(work)) then - ww => work(1:isz) - tx => work(isz+1:2*isz) - ty => work(2*isz+1:3*isz) - aux => work(3*isz+1:) - else if ((4*isz) <= size(work)) then - aux => work(1:) - allocate(ww(isz),tx(isz),ty(isz),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_alloc_request_,name,& - & i_err=(/3*isz,izero,izero,izero,izero/),& - & a_err='real(psb_dpk_)') - goto 9999 - end if - else if ((3*isz) <= size(work)) then - ww => work(1:isz) - tx => work(isz+1:2*isz) - ty => work(2*isz+1:3*isz) - allocate(aux(4*isz),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_alloc_request_,name,& - & i_err=(/4*isz,izero,izero,izero,izero/),& - & a_err='real(psb_dpk_)') - goto 9999 - end if + if (4*isz <= size(work)) then + aux => work(:) else - allocate(ww(isz),tx(isz),ty(isz),& - &aux(4*isz),stat=info) + allocate(aux(4*isz),stat=info) if (info /= psb_success_) then call psb_errpush(psb_err_alloc_request_,name,& & i_err=(/4*isz,izero,izero,izero,izero/),& & a_err='real(psb_dpk_)') goto 9999 end if - endif - if (sweeps == 0) then - + if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then ! - ! K^0 = I - ! zero sweeps of any smoother is just the identity. - ! - call psb_geaxpby(alpha,x,beta,y,desc_data,info) - - else if ((sm%novr == 0).and.(sweeps == 1).and.(.not.sm%sv%is_iterative())) then - ! - ! Shortcut: in this case it's just the same - ! as Block Jacobi. Moreover, if .not.sv%is_iterative, there's no need to pass init + ! Shortcut: in this case there is nothing else to be done. ! call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) @@ -153,388 +118,121 @@ subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& goto 9999 endif - else + else if (sweeps >= 0) then + ! + ! + ! Apply multiple sweeps of an AS solver + ! to compute an approximate solution of a linear system. + ! + ! + call psb_geasb(tx,sm%desc_data,info,mold=x%v,scratch=.true.) + call psb_geasb(ty,sm%desc_data,info,mold=x%v,scratch=.true.) + call psb_geasb(ww,sm%desc_data,info,mold=x%v,scratch=.true.) - call psb_geasb(vtx,sm%desc_data,info,mold=x%v,scratch=.true.) - call psb_geasb(vty,sm%desc_data,info,mold=x%v,scratch=.true.) - call psb_geasb(vww,sm%desc_data,info,mold=x%v,scratch=.true.) + ! + ! Unroll the first iteration and fold it inside SELECT CASE + ! this will save one SPMM when INIT=Z, and will be + ! significant when sweeps=1 (a common case) + ! + call psb_geaxpby(done,x,dzero,tx,desc_data,info) + if (info == 0) call sm%apply_restr(tx,trans_,aux,info) + if (info == 0) call psb_geaxpby(done,tx,dzero,ww,sm%desc_data,info) select case (init_) - case('Z') - call vtx%zero() + case('Z') + call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Z') + case('Y') - call psb_geaxpby(done,y,dzero,vtx,desc_data,info) + call psb_geaxpby(done,y,dzero,ty,desc_data,info) + if (info == 0) call sm%apply_restr(ty,trans_,aux,info) + if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,& + & work=aux,trans=trans_) + call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,info,init='Y') + case('U') if (.not.present(initu)) then call psb_errpush(psb_err_internal_error_,name,& & a_err='missing initu to smoother_apply') goto 9999 end if - call psb_geaxpby(done,initu,dzero,vtx,desc_data,info) + call psb_geaxpby(done,initu,dzero,ty,desc_data,info) + if (info == 0) call sm%apply_restr(ty,trans_,aux,info) + if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,& + & work=aux,trans=trans_) + call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,info,init='Y') + case default call psb_errpush(psb_err_internal_error_,name,& & a_err='wrong init to smoother_apply') goto 9999 end select - - if (sweeps == 1) then - - select case(trans_) - case('N') - ! - ! Get the overlap entries of tx (tx == x) - ! - if (sm%restr == psb_halo_) then - call psb_halo(vtx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - else if (sm%restr /= psb_none_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_restr_') - goto 9999 - end if - - - case('T','C') - ! - ! With transpose, we have to do it here - ! - - select case (sm%prol) - - case(psb_none_) - ! - ! Do nothing - - case(psb_sum_) - ! - ! The transpose of sum is halo - ! - call psb_halo(vtx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - - case(psb_avg_) - ! - ! Tricky one: first we have to scale the overlap entries, - ! which we can do by assignind mode=0, i.e. no communication - ! (hence only scaling), then we do the halo - ! - call psb_ovrl(vtx,sm%desc_data,info,& - & update=psb_avg_,work=aux,mode=izero) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - call psb_halo(vtx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_prol_') - goto 9999 - end select - - - case default - info=psb_err_iarg_invalid_i_ - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - - call sm%sv%apply(done,vtx,dzero,vty,sm%desc_data,trans_,aux,info,init='Y') - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in sub_aply Jacobi Sweeps = 1') - goto 9999 - endif - - select case(trans_) - case('N') - - select case (sm%prol) - - case(psb_none_) - ! - ! Would work anyway, but since it is supposed to do nothing ... - ! call psb_ovrl(ty,sm%desc_data,info,& - ! & update=sm%prol,work=aux) - - - case(psb_sum_,psb_avg_) - ! - ! Update the overlap of ty - ! - call psb_ovrl(vty,sm%desc_data,info,& - & update=sm%prol,work=aux) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_prol_') - goto 9999 - end select - - case('T','C') - ! - ! With transpose, we have to do it here - ! - if (sm%restr == psb_halo_) then - call psb_ovrl(vty,sm%desc_data,info,& - & update=psb_sum_,work=aux) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - else if (sm%restr /= psb_none_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_restr_') - goto 9999 - end if - - case default - info=psb_err_iarg_invalid_i_ - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - - - - else if (sweeps > 1) then - - ! - ! - ! Apply prec%iprcparm(mld_smoother_sweeps_) sweeps of a block-Jacobi solver - ! to compute an approximate solution of a linear system. + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in sub_aply Jacobi Sweeps = 1') + goto 9999 + endif + + do i=1, sweeps-1 ! + ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the + ! block diagonal part and the remaining part of the local matrix + ! and Y(j) is the approximate solution at sweep j. ! - select case (init_) - case('Z') - call vty%zero() - case('Y') - call psb_geaxpby(done,y,dzero,vty,sm%desc_data,info) - case('U') - if (.not.present(initu)) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='missing initu to smoother_apply') - goto 9999 - end if - call psb_geaxpby(done,initu,dzero,vty,sm%desc_data,info) - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='wrong init to smoother_apply') - goto 9999 - end select + if (info == 0) call psb_geaxpby(done,tx,dzero,ww,sm%desc_data,info) + if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,& + & work=aux,trans=trans_) - do i=1, sweeps - select case(trans_) - case('N') - ! - ! Get the overlap entries of tx (tx == x) - ! - if (sm%restr == psb_halo_) then - call psb_halo(vtx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - else if (sm%restr /= psb_none_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_restr_') - goto 9999 - end if - - - case('T','C') - ! - ! With transpose, we have to do it here - ! - - select case (sm%prol) - - case(psb_none_) - ! - ! Do nothing - - case(psb_sum_) - ! - ! The transpose of sum is halo - ! - call psb_halo(vtx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - - case(psb_avg_) - ! - ! Tricky one: first we have to scale the overlap entries, - ! which we can do by assignind mode=0, i.e. no communication - ! (hence only scaling), then we do the halo - ! - call psb_ovrl(vtx,sm%desc_data,info,& - & update=psb_avg_,work=aux,mode=izero) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - call psb_halo(vtx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_prol_') - goto 9999 - end select - - - case default - info=psb_err_iarg_invalid_i_ - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - ! - ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the - ! block diagonal part and the remaining part of the local matrix - ! and Y(j) is the approximate solution at sweep j. - ! - call psb_geaxpby(done,vtx,dzero,vww,sm%desc_data,info) - call psb_spmm(-done,sm%nd,vty,done,vww,sm%desc_data,info,& - & work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call sm%sv%apply(done,vww,dzero,vty,sm%desc_data,trans_,aux,info,init='Y') - - if (info /= psb_success_) exit - - - select case(trans_) - case('N') - - select case (sm%prol) - - case(psb_none_) - ! - ! Would work anyway, but since it is supposed to do nothing ... - ! call psb_ovrl(ty,sm%desc_data,info,& - ! & update=sm%prol,work=aux) - - - case(psb_sum_,psb_avg_) - ! - ! Update the overlap of ty - ! - call psb_ovrl(vty,sm%desc_data,info,& - & update=sm%prol,work=aux) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_prol_') - goto 9999 - end select + if (info /= psb_success_) exit - case('T','C') - ! - ! With transpose, we have to do it here - ! - if (sm%restr == psb_halo_) then - call psb_ovrl(vty,sm%desc_data,info,& - & update=psb_sum_,work=aux) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - else if (sm%restr /= psb_none_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_restr_') - goto 9999 - end if - - case default - info=psb_err_iarg_invalid_i_ - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - end do - - if (info /= psb_success_) then - info=psb_err_internal_error_ - call psb_errpush(info,name,& - & a_err='subsolve with Jacobi sweeps > 1') - goto 9999 - end if - - - else + call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Y') + + if (info /= psb_success_) exit + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + + end do - info = psb_err_iarg_neg_ + if (info /= psb_success_) then + info=psb_err_internal_error_ call psb_errpush(info,name,& - & i_err=(/itwo,sweeps,izero,izero,izero/)) - goto 9999 - - + & a_err='subsolve with Jacobi sweeps > 1') + goto 9999 end if - + ! ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) ! - call psb_geaxpby(alpha,vty,beta,y,desc_data,info) - - end if + call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + + + else + + info = psb_err_iarg_neg_ + call psb_errpush(info,name,& + & i_err=(/itwo,sweeps,izero,izero,izero/)) + goto 9999 + endif - if ((6*isz) <= size(work)) then - else if ((4*isz) <= size(work)) then - deallocate(ww,tx,ty) - else if ((3*isz) <= size(work)) then - deallocate(aux) - else - deallocate(ww,aux,tx,ty) + + if (.not.(4*isz <= size(work))) then + deallocate(aux,stat=info) endif - call vww%free(info) - call vtx%free(info) - call vty%free(info) + if (info ==0) call ww%free(info) + if (info ==0) call tx%free(info) + if (info ==0) call ty%free(info) + if (info /= 0) then + info = psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + end if call psb_erractionrestore(err_act) return - + 9999 call psb_error_handler(err_act) - + return end subroutine mld_d_as_smoother_apply_vect diff --git a/mlprec/impl/smoother/mld_d_as_smoother_prol_a.f90 b/mlprec/impl/smoother/mld_d_as_smoother_prol_a.f90 new file mode 100644 index 00000000..9b888bec --- /dev/null +++ b/mlprec/impl/smoother/mld_d_as_smoother_prol_a.f90 @@ -0,0 +1,150 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ 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_as_smoother_prol_a(sm,x,trans,work,info,data) + use psb_base_mod + use mld_d_as_smoother, mld_protect_nam => mld_d_as_smoother_prol_a + implicit none + class(mld_d_as_smoother_type), intent(inout) :: sm + real(psb_dpk_), intent(inout) :: x(:) + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: data + !Local + integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5), data_ + character :: trans_ + character(len=20) :: name='d_as_smther_prol_a', ch_err + + call psb_erractionsave(err_act) + + info = psb_success_ + ictxt = sm%desc_data%get_context() + call psb_info(ictxt,me,np) + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T') + case('C') + case default + info = psb_err_iarg_invalid_i_ + call psb_errpush(info,name) + goto 9999 + end select + + if (present(data)) then + data_ = data + else + data_ = psb_comm_ext_ + end if + + if (.not.allocated(sm%sv)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + + + select case(trans_) + case('N') + + select case (sm%prol) + + case(psb_none_) + ! + ! Would work anyway, but since it is supposed to do nothing ... + ! call psb_ovrl(x,sm%desc_data,info,& + ! & update=sm%prol,work=work) + + + case(psb_sum_,psb_avg_) + ! + ! Update the overlap of x + ! + call psb_ovrl(x,sm%desc_data,info,& + & update=sm%prol,work=work) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_ovrl' + goto 9999 + end if + + case default + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Invalid mld_sub_prol_') + goto 9999 + end select + + case('T','C') + ! + ! With transpose, we have to do it here + ! + if (sm%restr == psb_halo_) then + call psb_ovrl(x,sm%desc_data,info,& + & update=psb_sum_,work=work) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_ovrl' + goto 9999 + end if + else if (sm%restr /= psb_none_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Invalid mld_sub_restr_') + goto 9999 + end if + + case default + info=psb_err_iarg_invalid_i_ + int_err(1)=6 + ch_err(2:2)=trans + goto 9999 + end select + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_d_as_smoother_prol_a + + diff --git a/mlprec/impl/smoother/mld_d_as_smoother_prol_v.f90 b/mlprec/impl/smoother/mld_d_as_smoother_prol_v.f90 new file mode 100644 index 00000000..25ec0ff2 --- /dev/null +++ b/mlprec/impl/smoother/mld_d_as_smoother_prol_v.f90 @@ -0,0 +1,150 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ 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_as_smoother_prol_v(sm,x,trans,work,info,data) + use psb_base_mod + use mld_d_as_smoother, mld_protect_nam => mld_d_as_smoother_prol_v + implicit none + class(mld_d_as_smoother_type), intent(inout) :: sm + type(psb_d_vect_type),intent(inout) :: x + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: data + !Local + integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5), data_ + character :: trans_ + character(len=20) :: name='d_as_smther_prol_v', ch_err + + call psb_erractionsave(err_act) + + info = psb_success_ + ictxt = sm%desc_data%get_context() + call psb_info(ictxt,me,np) + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T') + case('C') + case default + info = psb_err_iarg_invalid_i_ + call psb_errpush(info,name) + goto 9999 + end select + + if (present(data)) then + data_ = data + else + data_ = psb_comm_ext_ + end if + + if (.not.allocated(sm%sv)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + + + select case(trans_) + case('N') + + select case (sm%prol) + + case(psb_none_) + ! + ! Would work anyway, but since it is supposed to do nothing ... + ! call psb_ovrl(x,sm%desc_data,info,& + ! & update=sm%prol,work=work) + + + case(psb_sum_,psb_avg_) + ! + ! Update the overlap of x + ! + call psb_ovrl(x,sm%desc_data,info,& + & update=sm%prol,work=work) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_ovrl' + goto 9999 + end if + + case default + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Invalid mld_sub_prol_') + goto 9999 + end select + + case('T','C') + ! + ! With transpose, we have to do it here + ! + if (sm%restr == psb_halo_) then + call psb_ovrl(x,sm%desc_data,info,& + & update=psb_sum_,work=work) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_ovrl' + goto 9999 + end if + else if (sm%restr /= psb_none_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Invalid mld_sub_restr_') + goto 9999 + end if + + case default + info=psb_err_iarg_invalid_i_ + int_err(1)=6 + ch_err(2:2)=trans + goto 9999 + end select + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_d_as_smoother_prol_v + + diff --git a/mlprec/impl/smoother/mld_d_as_smoother_restr_a.f90 b/mlprec/impl/smoother/mld_d_as_smoother_restr_a.f90 new file mode 100644 index 00000000..81d4fe3f --- /dev/null +++ b/mlprec/impl/smoother/mld_d_as_smoother_restr_a.f90 @@ -0,0 +1,169 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ 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_as_smoother_restr_a(sm,x,trans,work,info,data) + use psb_base_mod + use mld_d_as_smoother, mld_protect_nam => mld_d_as_smoother_restr_a + implicit none + class(mld_d_as_smoother_type), intent(inout) :: sm + real(psb_dpk_), intent(inout) :: x(:) + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: data + !Local + integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5), data_ + character :: trans_ + character(len=20) :: name='d_as_smther_restr_a', ch_err + + call psb_erractionsave(err_act) + + info = psb_success_ + ictxt = sm%desc_data%get_context() + call psb_info(ictxt,me,np) + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T') + case('C') + case default + info = psb_err_iarg_invalid_i_ + call psb_errpush(info,name) + goto 9999 + end select + + if (present(data)) then + data_ = data + else + data_ = psb_comm_ext_ + end if + + if (.not.allocated(sm%sv)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + + select case(trans_) + case('N') + ! + ! Get the overlap entries x + ! + if (sm%restr == psb_halo_) then + call psb_halo(x,sm%desc_data,info,work=work,data=data_) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_halo' + goto 9999 + end if + else if (sm%restr /= psb_none_) then + call psb_errpush(psb_err_internal_error_,name,& + &a_err='Invalid mld_sub_restr_') + goto 9999 + end if + + + case('T','C') + ! + ! With transpose, we have to do it here + ! + + select case (sm%prol) + + case(psb_none_) + ! + ! Do nothing + + case(psb_sum_) + ! + ! The transpose of sum is halo + ! + call psb_halo(x,sm%desc_data,info,work=work,data=data_) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_halo' + goto 9999 + end if + + case(psb_avg_) + ! + ! Tricky one: first we have to scale the overlap entries, + ! which we can do by assignind mode=0, i.e. no communication + ! (hence only scaling), then we do the halo + ! + call psb_ovrl(x,sm%desc_data,info,& + & update=psb_avg_,work=work,mode=izero) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_ovrl' + goto 9999 + end if + call psb_halo(x,sm%desc_data,info,work=work,data=data_) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_halo' + goto 9999 + end if + + case default + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Invalid mld_sub_prol_') + goto 9999 + end select + + + case default + info=psb_err_iarg_invalid_i_ + int_err(1)=6 + ch_err(2:2)=trans + goto 9999 + end select + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_d_as_smoother_restr_a + + diff --git a/mlprec/impl/smoother/mld_d_as_smoother_restr_v.f90 b/mlprec/impl/smoother/mld_d_as_smoother_restr_v.f90 new file mode 100644 index 00000000..02204ec7 --- /dev/null +++ b/mlprec/impl/smoother/mld_d_as_smoother_restr_v.f90 @@ -0,0 +1,169 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ 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_as_smoother_restr_v(sm,x,trans,work,info,data) + use psb_base_mod + use mld_d_as_smoother, mld_protect_nam => mld_d_as_smoother_restr_v + implicit none + class(mld_d_as_smoother_type), intent(inout) :: sm + type(psb_d_vect_type),intent(inout) :: x + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: data + !Local + integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5), data_ + character :: trans_ + character(len=20) :: name='d_as_smther_restr_v', ch_err + + call psb_erractionsave(err_act) + + info = psb_success_ + ictxt = sm%desc_data%get_context() + call psb_info(ictxt,me,np) + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T') + case('C') + case default + info = psb_err_iarg_invalid_i_ + call psb_errpush(info,name) + goto 9999 + end select + + if (present(data)) then + data_ = data + else + data_ = psb_comm_ext_ + end if + + if (.not.allocated(sm%sv)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + + select case(trans_) + case('N') + ! + ! Get the overlap entries x + ! + if (sm%restr == psb_halo_) then + call psb_halo(x,sm%desc_data,info,work=work,data=data_) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_halo' + goto 9999 + end if + else if (sm%restr /= psb_none_) then + call psb_errpush(psb_err_internal_error_,name,& + &a_err='Invalid mld_sub_restr_') + goto 9999 + end if + + + case('T','C') + ! + ! With transpose, we have to do it here + ! + + select case (sm%prol) + + case(psb_none_) + ! + ! Do nothing + + case(psb_sum_) + ! + ! The transpose of sum is halo + ! + call psb_halo(x,sm%desc_data,info,work=work,data=data_) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_halo' + goto 9999 + end if + + case(psb_avg_) + ! + ! Tricky one: first we have to scale the overlap entries, + ! which we can do by assignind mode=0, i.e. no communication + ! (hence only scaling), then we do the halo + ! + call psb_ovrl(x,sm%desc_data,info,& + & update=psb_avg_,work=work,mode=izero) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_ovrl' + goto 9999 + end if + call psb_halo(x,sm%desc_data,info,work=work,data=data_) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_halo' + goto 9999 + end if + + case default + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Invalid mld_sub_prol_') + goto 9999 + end select + + + case default + info=psb_err_iarg_invalid_i_ + int_err(1)=6 + ch_err(2:2)=trans + goto 9999 + end select + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_d_as_smoother_restr_v + + diff --git a/mlprec/impl/smoother/mld_s_as_smoother_apply.f90 b/mlprec/impl/smoother/mld_s_as_smoother_apply.f90 index 0033da63..410dac53 100644 --- a/mlprec/impl/smoother/mld_s_as_smoother_apply.f90 +++ b/mlprec/impl/smoother/mld_s_as_smoother_apply.f90 @@ -54,15 +54,15 @@ subroutine mld_s_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,& real(psb_spk_),intent(inout), optional :: initu(:) integer(psb_ipk_) :: n_row,n_col, nrow_d, i - real(psb_spk_), pointer :: ww(:), aux(:) - real(psb_spk_), allocatable :: tx(:),ty(:) + real(psb_spk_), pointer :: aux(:) + real(psb_spk_), allocatable :: tx(:),ty(:), ww(:) integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5) character :: trans_, init_ character(len=20) :: name='s_as_smoother_apply', ch_err call psb_erractionsave(err_act) - info = psb_success_ + info = psb_success_ ictxt = desc_data%get_context() call psb_info(ictxt,me,np) @@ -89,25 +89,14 @@ subroutine mld_s_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,& end if - n_row = sm%desc_data%get_local_rows() - n_col = sm%desc_data%get_local_cols() + n_row = sm%desc_data%get_local_rows() + n_col = sm%desc_data%get_local_cols() nrow_d = desc_data%get_local_rows() isz = max(n_row,N_COL) - if ((6*isz) <= size(work)) then - ww => work(1:isz) - aux => work(3*isz+1:) - else if ((4*isz) <= size(work)) then + if ((4*isz) <= size(work)) then aux => work(1:) - allocate(ww(isz),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_alloc_request_,name,& - & i_err=(/3*isz,izero,izero,izero,izero/),& - & a_err='real(psb_spk_)') - goto 9999 - end if - else if ((3*isz) <= size(work)) then - ww => work(1:isz) + else allocate(aux(4*isz),stat=info) if (info /= psb_success_) then call psb_errpush(psb_err_alloc_request_,name,& @@ -115,29 +104,11 @@ subroutine mld_s_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,& & a_err='real(psb_spk_)') goto 9999 end if - else - allocate(ww(isz), aux(4*isz),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_alloc_request_,name,& - & i_err=(/4*isz,izero,izero,izero,izero/),& - & a_err='real(psb_spk_)') - goto 9999 - end if - endif - if (sweeps == 0) then - + if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then ! - ! K^0 = I - ! zero sweeps of any smoother is just the identity. - ! - call psb_geaxpby(alpha,x,beta,y,desc_data,info) - - else if ((sm%novr == 0).and.(sweeps == 1).and.(.not.sm%sv%is_iterative())) then - ! - ! Shortcut: in this case it's just the same - ! as Block Jacobi. Moreover, if .not.sv%is_iterative, there's no need to pass init + ! Shortcut: in this case there is nothing else to be done. ! call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) @@ -147,380 +118,114 @@ subroutine mld_s_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,& goto 9999 endif - else - + else if (sweeps >= 0) then + ! + ! + ! Apply multiple sweeps of an AS solver + ! to compute an approximate solution of a linear system. + ! + ! + call psb_geasb(tx,sm%desc_data,info) + call psb_geasb(ty,sm%desc_data,info) + call psb_geasb(ww,sm%desc_data,info) - call psb_geasb(tx,desc_data,info) - call psb_geasb(ty,desc_data,info) + ! + ! Unroll the first iteration and fold it inside SELECT CASE + ! this will save one SPMM when INIT=Z, and will be + ! significant when sweeps=1 (a common case) + ! + call psb_geaxpby(sone,x,szero,tx,desc_data,info) + if (info == 0) call sm%apply_restr(tx,trans_,aux,info) + if (info == 0) call psb_geaxpby(sone,tx,szero,ww,sm%desc_data,info) select case (init_) - case('Z') - tx(:) = szero + case('Z') + call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Z') + case('Y') - call psb_geaxpby(sone,y,szero,tx,desc_data,info) + call psb_geaxpby(sone,y,szero,ty,desc_data,info) + if (info == 0) call sm%apply_restr(ty,trans_,aux,info) + if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,& + & work=aux,trans=trans_) + call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,info,init='Y') + case('U') if (.not.present(initu)) then call psb_errpush(psb_err_internal_error_,name,& & a_err='missing initu to smoother_apply') goto 9999 end if - call psb_geaxpby(sone,initu,szero,tx,desc_data,info) + call psb_geaxpby(sone,initu,szero,ty,desc_data,info) + if (info == 0) call sm%apply_restr(ty,trans_,aux,info) + if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,& + & work=aux,trans=trans_) + call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,info,init='Y') + case default call psb_errpush(psb_err_internal_error_,name,& & a_err='wrong init to smoother_apply') goto 9999 end select - - - if (sweeps == 1) then - - select case(trans_) - case('N') - ! - ! Get the overlap entries of tx (tx == x) - ! - if (sm%restr == psb_halo_) then - call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - else if (sm%restr /= psb_none_) then - call psb_errpush(psb_err_internal_error_,name,& - &a_err='Invalid mld_sub_restr_') - goto 9999 - end if - - - case('T','C') - ! - ! With transpose, we have to do it here - ! - - select case (sm%prol) - - case(psb_none_) - ! - ! Do nothing - - case(psb_sum_) - ! - ! The transpose of sum is halo - ! - call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - - case(psb_avg_) - ! - ! Tricky one: first we have to scale the overlap entries, - ! which we can do by assignind mode=0, i.e. no communication - ! (hence only scaling), then we do the halo - ! - call psb_ovrl(tx,sm%desc_data,info,& - & update=psb_avg_,work=aux,mode=izero) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_prol_') - goto 9999 - end select - - - case default - info=psb_err_iarg_invalid_i_ - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - - call sm%sv%apply(sone,tx,szero,ty,sm%desc_data,trans_,aux,info,init='Y') - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in sub_aply Jacobi Sweeps = 1') - goto 9999 - endif - - select case(trans_) - case('N') - - select case (sm%prol) - - case(psb_none_) - ! - ! Would work anyway, but since it is supposed to do nothing ... - ! call psb_ovrl(ty,sm%desc_data,info,& - ! & update=sm%prol,work=aux) - - - case(psb_sum_,psb_avg_) - ! - ! Update the overlap of ty - ! - call psb_ovrl(ty,sm%desc_data,info,& - & update=sm%prol,work=aux) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_prol_') - goto 9999 - end select - - case('T','C') - ! - ! With transpose, we have to do it here - ! - if (sm%restr == psb_halo_) then - call psb_ovrl(ty,sm%desc_data,info,& - & update=psb_sum_,work=aux) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - else if (sm%restr /= psb_none_) then - call psb_errpush(psb_err_internal_error_,name,& - &a_err='Invalid mld_sub_restr_') - goto 9999 - end if - - case default - info=psb_err_iarg_invalid_i_ - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - - - - else if (sweeps > 1) then - - ! - ! - ! Apply prec%iprcparm(mld_smoother_sweeps_) sweeps of a block-Jacobi solver - ! to compute an approximate solution of a linear system. + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in sub_aply Jacobi Sweeps = 1') + goto 9999 + endif + + do i=1, sweeps-1 ! + ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the + ! block diagonal part and the remaining part of the local matrix + ! and Y(j) is the approximate solution at sweep j. ! - select case (init_) - case('Z') - ty = szero - case('Y') - call psb_geaxpby(sone,y,szero,ty,sm%desc_data,info) - case('U') - if (.not.present(initu)) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='missing initu to smoother_apply') - goto 9999 - end if - call psb_geaxpby(sone,initu,szero,ty,sm%desc_data,info) - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='wrong init to smoother_apply') - goto 9999 - end select + if (info == 0) call psb_geaxpby(sone,tx,szero,ww,sm%desc_data,info) + if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,& + & work=aux,trans=trans_) - do i=1, sweeps - select case(trans_) - case('N') - ! - ! Get the overlap entries of tx (tx == x) - ! - if (sm%restr == psb_halo_) then - call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - else if (sm%restr /= psb_none_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_restr_') - goto 9999 - end if - - - case('T','C') - ! - ! With transpose, we have to do it here - ! - - select case (sm%prol) - - case(psb_none_) - ! - ! Do nothing - - case(psb_sum_) - ! - ! The transpose of sum is halo - ! - call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - - case(psb_avg_) - ! - ! Tricky one: first we have to scale the overlap entries, - ! which we can do by assignind mode=0, i.e. no communication - ! (hence only scaling), then we do the halo - ! - call psb_ovrl(tx,sm%desc_data,info,& - & update=psb_avg_,work=aux,mode=izero) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_prol_') - goto 9999 - end select - - - case default - info=psb_err_iarg_invalid_i_ - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - ! - ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the - ! block diagonal part and the remaining part of the local matrix - ! and Y(j) is the approximate solution at sweep j. - ! - ww(1:n_row) = tx(1:n_row) - call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,& - & work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Y') - - if (info /= psb_success_) exit - - - select case(trans_) - case('N') - - select case (sm%prol) - - case(psb_none_) - ! - ! Would work anyway, but since it is supposed to do nothing ... - ! call psb_ovrl(ty,sm%desc_data,info,& - ! & update=sm%prol,work=aux) - + if (info /= psb_success_) exit - case(psb_sum_,psb_avg_) - ! - ! Update the overlap of ty - ! - call psb_ovrl(ty,sm%desc_data,info,& - & update=sm%prol,work=aux) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_prol_') - goto 9999 - end select - - case('T','C') - ! - ! With transpose, we have to do it here - ! - if (sm%restr == psb_halo_) then - call psb_ovrl(ty,sm%desc_data,info,& - & update=psb_sum_,work=aux) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - else if (sm%restr /= psb_none_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_restr_') - goto 9999 - end if - - case default - info=psb_err_iarg_invalid_i_ - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - end do - - if (info /= psb_success_) then - info=psb_err_internal_error_ - call psb_errpush(info,name,a_err='subsolve with Jacobi sweeps > 1') - goto 9999 - end if - - - else + call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Y') + + if (info /= psb_success_) exit + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + + end do - info = psb_err_iarg_neg_ + if (info /= psb_success_) then + info=psb_err_internal_error_ call psb_errpush(info,name,& - & i_err=(/itwo,sweeps,izero,izero,izero/)) - goto 9999 - - + & a_err='subsolve with Jacobi sweeps > 1') + goto 9999 end if - + ! ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) ! call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + + + else + + info = psb_err_iarg_neg_ + call psb_errpush(info,name,& + & i_err=(/itwo,sweeps,izero,izero,izero/)) + goto 9999 - end if - - - if ((6*isz) <= size(work)) then - else if ((4*isz) <= size(work)) then - deallocate(ww,tx,ty) - else if ((3*isz) <= size(work)) then - deallocate(aux) - else - deallocate(ww,aux,tx,ty) endif + + if (.not.(4*isz <= size(work))) then + deallocate(aux,stat=info) + endif + if (info ==0) deallocate(ww,tx,ty,stat=info) + if (info /= 0) then + info = psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + end if + call psb_erractionrestore(err_act) return diff --git a/mlprec/impl/smoother/mld_s_as_smoother_apply_vect.f90 b/mlprec/impl/smoother/mld_s_as_smoother_apply_vect.f90 index 73e583ba..1dd04aef 100644 --- a/mlprec/impl/smoother/mld_s_as_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_s_as_smoother_apply_vect.f90 @@ -54,12 +54,11 @@ subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& type(psb_s_vect_type),intent(inout), optional :: initu integer(psb_ipk_) :: n_row,n_col, nrow_d, i - real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) - real(psb_spk_), allocatable :: vx(:) - type(psb_s_vect_type) :: vtx, vty, vww + real(psb_spk_), pointer :: aux(:) + type(psb_s_vect_type) :: tx, ty, ww integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5) character :: trans_, init_ - character(len=20) :: name='s_as_smoother_apply', ch_err + character(len=20) :: name='s_as_smoother_apply_v', ch_err call psb_erractionsave(err_act) @@ -95,55 +94,21 @@ subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& nrow_d = desc_data%get_local_rows() isz = max(n_row,N_COL) - if ((6*isz) <= size(work)) then - ww => work(1:isz) - tx => work(isz+1:2*isz) - ty => work(2*isz+1:3*isz) - aux => work(3*isz+1:) - else if ((4*isz) <= size(work)) then - aux => work(1:) - allocate(ww(isz),tx(isz),ty(isz),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_alloc_request_,name,& - & i_err=(/3*isz,izero,izero,izero,izero/),& - & a_err='real(psb_spk_)') - goto 9999 - end if - else if ((3*isz) <= size(work)) then - ww => work(1:isz) - tx => work(isz+1:2*isz) - ty => work(2*isz+1:3*isz) - allocate(aux(4*isz),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_alloc_request_,name,& - & i_err=(/4*isz,izero,izero,izero,izero/),& - & a_err='real(psb_spk_)') - goto 9999 - end if + if (4*isz <= size(work)) then + aux => work(:) else - allocate(ww(isz),tx(isz),ty(isz),& - &aux(4*isz),stat=info) + allocate(aux(4*isz),stat=info) if (info /= psb_success_) then call psb_errpush(psb_err_alloc_request_,name,& & i_err=(/4*isz,izero,izero,izero,izero/),& & a_err='real(psb_spk_)') goto 9999 end if - endif - if (sweeps == 0) then - + if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then ! - ! K^0 = I - ! zero sweeps of any smoother is just the identity. - ! - call psb_geaxpby(alpha,x,beta,y,desc_data,info) - - else if ((sm%novr == 0).and.(sweeps == 1).and.(.not.sm%sv%is_iterative())) then - ! - ! Shortcut: in this case it's just the same - ! as Block Jacobi. Moreover, if .not.sv%is_iterative, there's no need to pass init + ! Shortcut: in this case there is nothing else to be done. ! call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) @@ -153,388 +118,121 @@ subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& goto 9999 endif - else + else if (sweeps >= 0) then + ! + ! + ! Apply multiple sweeps of an AS solver + ! to compute an approximate solution of a linear system. + ! + ! + call psb_geasb(tx,sm%desc_data,info,mold=x%v,scratch=.true.) + call psb_geasb(ty,sm%desc_data,info,mold=x%v,scratch=.true.) + call psb_geasb(ww,sm%desc_data,info,mold=x%v,scratch=.true.) - call psb_geasb(vtx,sm%desc_data,info,mold=x%v,scratch=.true.) - call psb_geasb(vty,sm%desc_data,info,mold=x%v,scratch=.true.) - call psb_geasb(vww,sm%desc_data,info,mold=x%v,scratch=.true.) + ! + ! Unroll the first iteration and fold it inside SELECT CASE + ! this will save one SPMM when INIT=Z, and will be + ! significant when sweeps=1 (a common case) + ! + call psb_geaxpby(sone,x,szero,tx,desc_data,info) + if (info == 0) call sm%apply_restr(tx,trans_,aux,info) + if (info == 0) call psb_geaxpby(sone,tx,szero,ww,sm%desc_data,info) select case (init_) - case('Z') - call vtx%zero() + case('Z') + call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Z') + case('Y') - call psb_geaxpby(sone,y,szero,vtx,desc_data,info) + call psb_geaxpby(sone,y,szero,ty,desc_data,info) + if (info == 0) call sm%apply_restr(ty,trans_,aux,info) + if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,& + & work=aux,trans=trans_) + call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,info,init='Y') + case('U') if (.not.present(initu)) then call psb_errpush(psb_err_internal_error_,name,& & a_err='missing initu to smoother_apply') goto 9999 end if - call psb_geaxpby(sone,initu,szero,vtx,desc_data,info) + call psb_geaxpby(sone,initu,szero,ty,desc_data,info) + if (info == 0) call sm%apply_restr(ty,trans_,aux,info) + if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,& + & work=aux,trans=trans_) + call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,info,init='Y') + case default call psb_errpush(psb_err_internal_error_,name,& & a_err='wrong init to smoother_apply') goto 9999 end select - - if (sweeps == 1) then - - select case(trans_) - case('N') - ! - ! Get the overlap entries of tx (tx == x) - ! - if (sm%restr == psb_halo_) then - call psb_halo(vtx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - else if (sm%restr /= psb_none_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_restr_') - goto 9999 - end if - - - case('T','C') - ! - ! With transpose, we have to do it here - ! - - select case (sm%prol) - - case(psb_none_) - ! - ! Do nothing - - case(psb_sum_) - ! - ! The transpose of sum is halo - ! - call psb_halo(vtx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - - case(psb_avg_) - ! - ! Tricky one: first we have to scale the overlap entries, - ! which we can do by assignind mode=0, i.e. no communication - ! (hence only scaling), then we do the halo - ! - call psb_ovrl(vtx,sm%desc_data,info,& - & update=psb_avg_,work=aux,mode=izero) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - call psb_halo(vtx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_prol_') - goto 9999 - end select - - - case default - info=psb_err_iarg_invalid_i_ - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - - call sm%sv%apply(sone,vtx,szero,vty,sm%desc_data,trans_,aux,info,init='Y') - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in sub_aply Jacobi Sweeps = 1') - goto 9999 - endif - - select case(trans_) - case('N') - - select case (sm%prol) - - case(psb_none_) - ! - ! Would work anyway, but since it is supposed to do nothing ... - ! call psb_ovrl(ty,sm%desc_data,info,& - ! & update=sm%prol,work=aux) - - - case(psb_sum_,psb_avg_) - ! - ! Update the overlap of ty - ! - call psb_ovrl(vty,sm%desc_data,info,& - & update=sm%prol,work=aux) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_prol_') - goto 9999 - end select - - case('T','C') - ! - ! With transpose, we have to do it here - ! - if (sm%restr == psb_halo_) then - call psb_ovrl(vty,sm%desc_data,info,& - & update=psb_sum_,work=aux) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - else if (sm%restr /= psb_none_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_restr_') - goto 9999 - end if - - case default - info=psb_err_iarg_invalid_i_ - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - - - - else if (sweeps > 1) then - - ! - ! - ! Apply prec%iprcparm(mld_smoother_sweeps_) sweeps of a block-Jacobi solver - ! to compute an approximate solution of a linear system. + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in sub_aply Jacobi Sweeps = 1') + goto 9999 + endif + + do i=1, sweeps-1 ! + ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the + ! block diagonal part and the remaining part of the local matrix + ! and Y(j) is the approximate solution at sweep j. ! - select case (init_) - case('Z') - call vty%zero() - case('Y') - call psb_geaxpby(sone,y,szero,vty,sm%desc_data,info) - case('U') - if (.not.present(initu)) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='missing initu to smoother_apply') - goto 9999 - end if - call psb_geaxpby(sone,initu,szero,vty,sm%desc_data,info) - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='wrong init to smoother_apply') - goto 9999 - end select + if (info == 0) call psb_geaxpby(sone,tx,szero,ww,sm%desc_data,info) + if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,& + & work=aux,trans=trans_) - do i=1, sweeps - select case(trans_) - case('N') - ! - ! Get the overlap entries of tx (tx == x) - ! - if (sm%restr == psb_halo_) then - call psb_halo(vtx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - else if (sm%restr /= psb_none_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_restr_') - goto 9999 - end if - - - case('T','C') - ! - ! With transpose, we have to do it here - ! - - select case (sm%prol) - - case(psb_none_) - ! - ! Do nothing - - case(psb_sum_) - ! - ! The transpose of sum is halo - ! - call psb_halo(vtx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - - case(psb_avg_) - ! - ! Tricky one: first we have to scale the overlap entries, - ! which we can do by assignind mode=0, i.e. no communication - ! (hence only scaling), then we do the halo - ! - call psb_ovrl(vtx,sm%desc_data,info,& - & update=psb_avg_,work=aux,mode=izero) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - call psb_halo(vtx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_prol_') - goto 9999 - end select - - - case default - info=psb_err_iarg_invalid_i_ - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - ! - ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the - ! block diagonal part and the remaining part of the local matrix - ! and Y(j) is the approximate solution at sweep j. - ! - call psb_geaxpby(sone,vtx,szero,vww,sm%desc_data,info) - call psb_spmm(-sone,sm%nd,vty,sone,vww,sm%desc_data,info,& - & work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call sm%sv%apply(sone,vww,szero,vty,sm%desc_data,trans_,aux,info,init='Y') - - if (info /= psb_success_) exit - - - select case(trans_) - case('N') - - select case (sm%prol) - - case(psb_none_) - ! - ! Would work anyway, but since it is supposed to do nothing ... - ! call psb_ovrl(ty,sm%desc_data,info,& - ! & update=sm%prol,work=aux) - - - case(psb_sum_,psb_avg_) - ! - ! Update the overlap of ty - ! - call psb_ovrl(vty,sm%desc_data,info,& - & update=sm%prol,work=aux) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_prol_') - goto 9999 - end select + if (info /= psb_success_) exit - case('T','C') - ! - ! With transpose, we have to do it here - ! - if (sm%restr == psb_halo_) then - call psb_ovrl(vty,sm%desc_data,info,& - & update=psb_sum_,work=aux) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - else if (sm%restr /= psb_none_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_restr_') - goto 9999 - end if - - case default - info=psb_err_iarg_invalid_i_ - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - end do - - if (info /= psb_success_) then - info=psb_err_internal_error_ - call psb_errpush(info,name,& - & a_err='subsolve with Jacobi sweeps > 1') - goto 9999 - end if - - - else + call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Y') + + if (info /= psb_success_) exit + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + + end do - info = psb_err_iarg_neg_ + if (info /= psb_success_) then + info=psb_err_internal_error_ call psb_errpush(info,name,& - & i_err=(/itwo,sweeps,izero,izero,izero/)) - goto 9999 - - + & a_err='subsolve with Jacobi sweeps > 1') + goto 9999 end if - + ! ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) ! - call psb_geaxpby(alpha,vty,beta,y,desc_data,info) - - end if + call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + + + else + + info = psb_err_iarg_neg_ + call psb_errpush(info,name,& + & i_err=(/itwo,sweeps,izero,izero,izero/)) + goto 9999 + endif - if ((6*isz) <= size(work)) then - else if ((4*isz) <= size(work)) then - deallocate(ww,tx,ty) - else if ((3*isz) <= size(work)) then - deallocate(aux) - else - deallocate(ww,aux,tx,ty) + + if (.not.(4*isz <= size(work))) then + deallocate(aux,stat=info) endif - call vww%free(info) - call vtx%free(info) - call vty%free(info) + if (info ==0) call ww%free(info) + if (info ==0) call tx%free(info) + if (info ==0) call ty%free(info) + if (info /= 0) then + info = psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + end if call psb_erractionrestore(err_act) return - + 9999 call psb_error_handler(err_act) - + return end subroutine mld_s_as_smoother_apply_vect diff --git a/mlprec/impl/smoother/mld_s_as_smoother_prol_a.f90 b/mlprec/impl/smoother/mld_s_as_smoother_prol_a.f90 new file mode 100644 index 00000000..56a38878 --- /dev/null +++ b/mlprec/impl/smoother/mld_s_as_smoother_prol_a.f90 @@ -0,0 +1,150 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ 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_as_smoother_prol_a(sm,x,trans,work,info,data) + use psb_base_mod + use mld_s_as_smoother, mld_protect_nam => mld_s_as_smoother_prol_a + implicit none + class(mld_s_as_smoother_type), intent(inout) :: sm + real(psb_spk_), intent(inout) :: x(:) + character(len=1),intent(in) :: trans + real(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: data + !Local + integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5), data_ + character :: trans_ + character(len=20) :: name='s_as_smther_prol_a', ch_err + + call psb_erractionsave(err_act) + + info = psb_success_ + ictxt = sm%desc_data%get_context() + call psb_info(ictxt,me,np) + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T') + case('C') + case default + info = psb_err_iarg_invalid_i_ + call psb_errpush(info,name) + goto 9999 + end select + + if (present(data)) then + data_ = data + else + data_ = psb_comm_ext_ + end if + + if (.not.allocated(sm%sv)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + + + select case(trans_) + case('N') + + select case (sm%prol) + + case(psb_none_) + ! + ! Would work anyway, but since it is supposed to do nothing ... + ! call psb_ovrl(x,sm%desc_data,info,& + ! & update=sm%prol,work=work) + + + case(psb_sum_,psb_avg_) + ! + ! Update the overlap of x + ! + call psb_ovrl(x,sm%desc_data,info,& + & update=sm%prol,work=work) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_ovrl' + goto 9999 + end if + + case default + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Invalid mld_sub_prol_') + goto 9999 + end select + + case('T','C') + ! + ! With transpose, we have to do it here + ! + if (sm%restr == psb_halo_) then + call psb_ovrl(x,sm%desc_data,info,& + & update=psb_sum_,work=work) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_ovrl' + goto 9999 + end if + else if (sm%restr /= psb_none_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Invalid mld_sub_restr_') + goto 9999 + end if + + case default + info=psb_err_iarg_invalid_i_ + int_err(1)=6 + ch_err(2:2)=trans + goto 9999 + end select + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_s_as_smoother_prol_a + + diff --git a/mlprec/impl/smoother/mld_s_as_smoother_prol_v.f90 b/mlprec/impl/smoother/mld_s_as_smoother_prol_v.f90 new file mode 100644 index 00000000..2b82d355 --- /dev/null +++ b/mlprec/impl/smoother/mld_s_as_smoother_prol_v.f90 @@ -0,0 +1,150 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ 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_as_smoother_prol_v(sm,x,trans,work,info,data) + use psb_base_mod + use mld_s_as_smoother, mld_protect_nam => mld_s_as_smoother_prol_v + implicit none + class(mld_s_as_smoother_type), intent(inout) :: sm + type(psb_s_vect_type),intent(inout) :: x + character(len=1),intent(in) :: trans + real(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: data + !Local + integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5), data_ + character :: trans_ + character(len=20) :: name='s_as_smther_prol_v', ch_err + + call psb_erractionsave(err_act) + + info = psb_success_ + ictxt = sm%desc_data%get_context() + call psb_info(ictxt,me,np) + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T') + case('C') + case default + info = psb_err_iarg_invalid_i_ + call psb_errpush(info,name) + goto 9999 + end select + + if (present(data)) then + data_ = data + else + data_ = psb_comm_ext_ + end if + + if (.not.allocated(sm%sv)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + + + select case(trans_) + case('N') + + select case (sm%prol) + + case(psb_none_) + ! + ! Would work anyway, but since it is supposed to do nothing ... + ! call psb_ovrl(x,sm%desc_data,info,& + ! & update=sm%prol,work=work) + + + case(psb_sum_,psb_avg_) + ! + ! Update the overlap of x + ! + call psb_ovrl(x,sm%desc_data,info,& + & update=sm%prol,work=work) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_ovrl' + goto 9999 + end if + + case default + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Invalid mld_sub_prol_') + goto 9999 + end select + + case('T','C') + ! + ! With transpose, we have to do it here + ! + if (sm%restr == psb_halo_) then + call psb_ovrl(x,sm%desc_data,info,& + & update=psb_sum_,work=work) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_ovrl' + goto 9999 + end if + else if (sm%restr /= psb_none_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Invalid mld_sub_restr_') + goto 9999 + end if + + case default + info=psb_err_iarg_invalid_i_ + int_err(1)=6 + ch_err(2:2)=trans + goto 9999 + end select + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_s_as_smoother_prol_v + + diff --git a/mlprec/impl/smoother/mld_s_as_smoother_restr_a.f90 b/mlprec/impl/smoother/mld_s_as_smoother_restr_a.f90 new file mode 100644 index 00000000..4337f74a --- /dev/null +++ b/mlprec/impl/smoother/mld_s_as_smoother_restr_a.f90 @@ -0,0 +1,169 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ 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_as_smoother_restr_a(sm,x,trans,work,info,data) + use psb_base_mod + use mld_s_as_smoother, mld_protect_nam => mld_s_as_smoother_restr_a + implicit none + class(mld_s_as_smoother_type), intent(inout) :: sm + real(psb_spk_), intent(inout) :: x(:) + character(len=1),intent(in) :: trans + real(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: data + !Local + integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5), data_ + character :: trans_ + character(len=20) :: name='s_as_smther_restr_a', ch_err + + call psb_erractionsave(err_act) + + info = psb_success_ + ictxt = sm%desc_data%get_context() + call psb_info(ictxt,me,np) + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T') + case('C') + case default + info = psb_err_iarg_invalid_i_ + call psb_errpush(info,name) + goto 9999 + end select + + if (present(data)) then + data_ = data + else + data_ = psb_comm_ext_ + end if + + if (.not.allocated(sm%sv)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + + select case(trans_) + case('N') + ! + ! Get the overlap entries x + ! + if (sm%restr == psb_halo_) then + call psb_halo(x,sm%desc_data,info,work=work,data=data_) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_halo' + goto 9999 + end if + else if (sm%restr /= psb_none_) then + call psb_errpush(psb_err_internal_error_,name,& + &a_err='Invalid mld_sub_restr_') + goto 9999 + end if + + + case('T','C') + ! + ! With transpose, we have to do it here + ! + + select case (sm%prol) + + case(psb_none_) + ! + ! Do nothing + + case(psb_sum_) + ! + ! The transpose of sum is halo + ! + call psb_halo(x,sm%desc_data,info,work=work,data=data_) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_halo' + goto 9999 + end if + + case(psb_avg_) + ! + ! Tricky one: first we have to scale the overlap entries, + ! which we can do by assignind mode=0, i.e. no communication + ! (hence only scaling), then we do the halo + ! + call psb_ovrl(x,sm%desc_data,info,& + & update=psb_avg_,work=work,mode=izero) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_ovrl' + goto 9999 + end if + call psb_halo(x,sm%desc_data,info,work=work,data=data_) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_halo' + goto 9999 + end if + + case default + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Invalid mld_sub_prol_') + goto 9999 + end select + + + case default + info=psb_err_iarg_invalid_i_ + int_err(1)=6 + ch_err(2:2)=trans + goto 9999 + end select + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_s_as_smoother_restr_a + + diff --git a/mlprec/impl/smoother/mld_s_as_smoother_restr_v.f90 b/mlprec/impl/smoother/mld_s_as_smoother_restr_v.f90 new file mode 100644 index 00000000..adaf7a1e --- /dev/null +++ b/mlprec/impl/smoother/mld_s_as_smoother_restr_v.f90 @@ -0,0 +1,169 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ 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_as_smoother_restr_v(sm,x,trans,work,info,data) + use psb_base_mod + use mld_s_as_smoother, mld_protect_nam => mld_s_as_smoother_restr_v + implicit none + class(mld_s_as_smoother_type), intent(inout) :: sm + type(psb_s_vect_type),intent(inout) :: x + character(len=1),intent(in) :: trans + real(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: data + !Local + integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5), data_ + character :: trans_ + character(len=20) :: name='s_as_smther_restr_v', ch_err + + call psb_erractionsave(err_act) + + info = psb_success_ + ictxt = sm%desc_data%get_context() + call psb_info(ictxt,me,np) + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T') + case('C') + case default + info = psb_err_iarg_invalid_i_ + call psb_errpush(info,name) + goto 9999 + end select + + if (present(data)) then + data_ = data + else + data_ = psb_comm_ext_ + end if + + if (.not.allocated(sm%sv)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + + select case(trans_) + case('N') + ! + ! Get the overlap entries x + ! + if (sm%restr == psb_halo_) then + call psb_halo(x,sm%desc_data,info,work=work,data=data_) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_halo' + goto 9999 + end if + else if (sm%restr /= psb_none_) then + call psb_errpush(psb_err_internal_error_,name,& + &a_err='Invalid mld_sub_restr_') + goto 9999 + end if + + + case('T','C') + ! + ! With transpose, we have to do it here + ! + + select case (sm%prol) + + case(psb_none_) + ! + ! Do nothing + + case(psb_sum_) + ! + ! The transpose of sum is halo + ! + call psb_halo(x,sm%desc_data,info,work=work,data=data_) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_halo' + goto 9999 + end if + + case(psb_avg_) + ! + ! Tricky one: first we have to scale the overlap entries, + ! which we can do by assignind mode=0, i.e. no communication + ! (hence only scaling), then we do the halo + ! + call psb_ovrl(x,sm%desc_data,info,& + & update=psb_avg_,work=work,mode=izero) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_ovrl' + goto 9999 + end if + call psb_halo(x,sm%desc_data,info,work=work,data=data_) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_halo' + goto 9999 + end if + + case default + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Invalid mld_sub_prol_') + goto 9999 + end select + + + case default + info=psb_err_iarg_invalid_i_ + int_err(1)=6 + ch_err(2:2)=trans + goto 9999 + end select + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_s_as_smoother_restr_v + + diff --git a/mlprec/impl/smoother/mld_z_as_smoother_apply.f90 b/mlprec/impl/smoother/mld_z_as_smoother_apply.f90 index b74aebba..17b2ae8a 100644 --- a/mlprec/impl/smoother/mld_z_as_smoother_apply.f90 +++ b/mlprec/impl/smoother/mld_z_as_smoother_apply.f90 @@ -54,15 +54,15 @@ subroutine mld_z_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,& complex(psb_dpk_),intent(inout), optional :: initu(:) integer(psb_ipk_) :: n_row,n_col, nrow_d, i - complex(psb_dpk_), pointer :: ww(:), aux(:) - complex(psb_dpk_), allocatable :: tx(:),ty(:) + complex(psb_dpk_), pointer :: aux(:) + complex(psb_dpk_), allocatable :: tx(:),ty(:), ww(:) integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5) character :: trans_, init_ character(len=20) :: name='z_as_smoother_apply', ch_err call psb_erractionsave(err_act) - info = psb_success_ + info = psb_success_ ictxt = desc_data%get_context() call psb_info(ictxt,me,np) @@ -89,25 +89,14 @@ subroutine mld_z_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,& end if - n_row = sm%desc_data%get_local_rows() - n_col = sm%desc_data%get_local_cols() + n_row = sm%desc_data%get_local_rows() + n_col = sm%desc_data%get_local_cols() nrow_d = desc_data%get_local_rows() isz = max(n_row,N_COL) - if ((6*isz) <= size(work)) then - ww => work(1:isz) - aux => work(3*isz+1:) - else if ((4*isz) <= size(work)) then + if ((4*isz) <= size(work)) then aux => work(1:) - allocate(ww(isz),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_alloc_request_,name,& - & i_err=(/3*isz,izero,izero,izero,izero/),& - & a_err='complex(psb_dpk_)') - goto 9999 - end if - else if ((3*isz) <= size(work)) then - ww => work(1:isz) + else allocate(aux(4*isz),stat=info) if (info /= psb_success_) then call psb_errpush(psb_err_alloc_request_,name,& @@ -115,29 +104,11 @@ subroutine mld_z_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,& & a_err='complex(psb_dpk_)') goto 9999 end if - else - allocate(ww(isz), aux(4*isz),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_alloc_request_,name,& - & i_err=(/4*isz,izero,izero,izero,izero/),& - & a_err='complex(psb_dpk_)') - goto 9999 - end if - endif - if (sweeps == 0) then - + if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then ! - ! K^0 = I - ! zero sweeps of any smoother is just the identity. - ! - call psb_geaxpby(alpha,x,beta,y,desc_data,info) - - else if ((sm%novr == 0).and.(sweeps == 1).and.(.not.sm%sv%is_iterative())) then - ! - ! Shortcut: in this case it's just the same - ! as Block Jacobi. Moreover, if .not.sv%is_iterative, there's no need to pass init + ! Shortcut: in this case there is nothing else to be done. ! call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) @@ -147,380 +118,114 @@ subroutine mld_z_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,& goto 9999 endif - else - + else if (sweeps >= 0) then + ! + ! + ! Apply multiple sweeps of an AS solver + ! to compute an approximate solution of a linear system. + ! + ! + call psb_geasb(tx,sm%desc_data,info) + call psb_geasb(ty,sm%desc_data,info) + call psb_geasb(ww,sm%desc_data,info) - call psb_geasb(tx,desc_data,info) - call psb_geasb(ty,desc_data,info) + ! + ! Unroll the first iteration and fold it inside SELECT CASE + ! this will save one SPMM when INIT=Z, and will be + ! significant when sweeps=1 (a common case) + ! + call psb_geaxpby(zone,x,zzero,tx,desc_data,info) + if (info == 0) call sm%apply_restr(tx,trans_,aux,info) + if (info == 0) call psb_geaxpby(zone,tx,zzero,ww,sm%desc_data,info) select case (init_) - case('Z') - tx(:) = zzero + case('Z') + call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Z') + case('Y') - call psb_geaxpby(zone,y,zzero,tx,desc_data,info) + call psb_geaxpby(zone,y,zzero,ty,desc_data,info) + if (info == 0) call sm%apply_restr(ty,trans_,aux,info) + if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,& + & work=aux,trans=trans_) + call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,info,init='Y') + case('U') if (.not.present(initu)) then call psb_errpush(psb_err_internal_error_,name,& & a_err='missing initu to smoother_apply') goto 9999 end if - call psb_geaxpby(zone,initu,zzero,tx,desc_data,info) + call psb_geaxpby(zone,initu,zzero,ty,desc_data,info) + if (info == 0) call sm%apply_restr(ty,trans_,aux,info) + if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,& + & work=aux,trans=trans_) + call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,info,init='Y') + case default call psb_errpush(psb_err_internal_error_,name,& & a_err='wrong init to smoother_apply') goto 9999 end select - - - if (sweeps == 1) then - - select case(trans_) - case('N') - ! - ! Get the overlap entries of tx (tx == x) - ! - if (sm%restr == psb_halo_) then - call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - else if (sm%restr /= psb_none_) then - call psb_errpush(psb_err_internal_error_,name,& - &a_err='Invalid mld_sub_restr_') - goto 9999 - end if - - - case('T','C') - ! - ! With transpose, we have to do it here - ! - - select case (sm%prol) - - case(psb_none_) - ! - ! Do nothing - - case(psb_sum_) - ! - ! The transpose of sum is halo - ! - call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - - case(psb_avg_) - ! - ! Tricky one: first we have to scale the overlap entries, - ! which we can do by assignind mode=0, i.e. no communication - ! (hence only scaling), then we do the halo - ! - call psb_ovrl(tx,sm%desc_data,info,& - & update=psb_avg_,work=aux,mode=izero) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_prol_') - goto 9999 - end select - - - case default - info=psb_err_iarg_invalid_i_ - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - - call sm%sv%apply(zone,tx,zzero,ty,sm%desc_data,trans_,aux,info,init='Y') - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in sub_aply Jacobi Sweeps = 1') - goto 9999 - endif - - select case(trans_) - case('N') - - select case (sm%prol) - - case(psb_none_) - ! - ! Would work anyway, but since it is supposed to do nothing ... - ! call psb_ovrl(ty,sm%desc_data,info,& - ! & update=sm%prol,work=aux) - - - case(psb_sum_,psb_avg_) - ! - ! Update the overlap of ty - ! - call psb_ovrl(ty,sm%desc_data,info,& - & update=sm%prol,work=aux) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_prol_') - goto 9999 - end select - - case('T','C') - ! - ! With transpose, we have to do it here - ! - if (sm%restr == psb_halo_) then - call psb_ovrl(ty,sm%desc_data,info,& - & update=psb_sum_,work=aux) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - else if (sm%restr /= psb_none_) then - call psb_errpush(psb_err_internal_error_,name,& - &a_err='Invalid mld_sub_restr_') - goto 9999 - end if - - case default - info=psb_err_iarg_invalid_i_ - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - - - - else if (sweeps > 1) then - - ! - ! - ! Apply prec%iprcparm(mld_smoother_sweeps_) sweeps of a block-Jacobi solver - ! to compute an approximate solution of a linear system. + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in sub_aply Jacobi Sweeps = 1') + goto 9999 + endif + + do i=1, sweeps-1 ! + ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the + ! block diagonal part and the remaining part of the local matrix + ! and Y(j) is the approximate solution at sweep j. ! - select case (init_) - case('Z') - ty = zzero - case('Y') - call psb_geaxpby(zone,y,zzero,ty,sm%desc_data,info) - case('U') - if (.not.present(initu)) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='missing initu to smoother_apply') - goto 9999 - end if - call psb_geaxpby(zone,initu,zzero,ty,sm%desc_data,info) - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='wrong init to smoother_apply') - goto 9999 - end select + if (info == 0) call psb_geaxpby(zone,tx,zzero,ww,sm%desc_data,info) + if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,& + & work=aux,trans=trans_) - do i=1, sweeps - select case(trans_) - case('N') - ! - ! Get the overlap entries of tx (tx == x) - ! - if (sm%restr == psb_halo_) then - call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - else if (sm%restr /= psb_none_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_restr_') - goto 9999 - end if - - - case('T','C') - ! - ! With transpose, we have to do it here - ! - - select case (sm%prol) - - case(psb_none_) - ! - ! Do nothing - - case(psb_sum_) - ! - ! The transpose of sum is halo - ! - call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - - case(psb_avg_) - ! - ! Tricky one: first we have to scale the overlap entries, - ! which we can do by assignind mode=0, i.e. no communication - ! (hence only scaling), then we do the halo - ! - call psb_ovrl(tx,sm%desc_data,info,& - & update=psb_avg_,work=aux,mode=izero) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_prol_') - goto 9999 - end select - - - case default - info=psb_err_iarg_invalid_i_ - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - ! - ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the - ! block diagonal part and the remaining part of the local matrix - ! and Y(j) is the approximate solution at sweep j. - ! - ww(1:n_row) = tx(1:n_row) - call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,& - & work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Y') - - if (info /= psb_success_) exit - - - select case(trans_) - case('N') - - select case (sm%prol) - - case(psb_none_) - ! - ! Would work anyway, but since it is supposed to do nothing ... - ! call psb_ovrl(ty,sm%desc_data,info,& - ! & update=sm%prol,work=aux) - + if (info /= psb_success_) exit - case(psb_sum_,psb_avg_) - ! - ! Update the overlap of ty - ! - call psb_ovrl(ty,sm%desc_data,info,& - & update=sm%prol,work=aux) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_prol_') - goto 9999 - end select - - case('T','C') - ! - ! With transpose, we have to do it here - ! - if (sm%restr == psb_halo_) then - call psb_ovrl(ty,sm%desc_data,info,& - & update=psb_sum_,work=aux) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - else if (sm%restr /= psb_none_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_restr_') - goto 9999 - end if - - case default - info=psb_err_iarg_invalid_i_ - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - end do - - if (info /= psb_success_) then - info=psb_err_internal_error_ - call psb_errpush(info,name,a_err='subsolve with Jacobi sweeps > 1') - goto 9999 - end if - - - else + call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Y') + + if (info /= psb_success_) exit + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + + end do - info = psb_err_iarg_neg_ + if (info /= psb_success_) then + info=psb_err_internal_error_ call psb_errpush(info,name,& - & i_err=(/itwo,sweeps,izero,izero,izero/)) - goto 9999 - - + & a_err='subsolve with Jacobi sweeps > 1') + goto 9999 end if - + ! ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) ! call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + + + else + + info = psb_err_iarg_neg_ + call psb_errpush(info,name,& + & i_err=(/itwo,sweeps,izero,izero,izero/)) + goto 9999 - end if - - - if ((6*isz) <= size(work)) then - else if ((4*isz) <= size(work)) then - deallocate(ww,tx,ty) - else if ((3*isz) <= size(work)) then - deallocate(aux) - else - deallocate(ww,aux,tx,ty) endif + + if (.not.(4*isz <= size(work))) then + deallocate(aux,stat=info) + endif + if (info ==0) deallocate(ww,tx,ty,stat=info) + if (info /= 0) then + info = psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + end if + call psb_erractionrestore(err_act) return diff --git a/mlprec/impl/smoother/mld_z_as_smoother_apply_vect.f90 b/mlprec/impl/smoother/mld_z_as_smoother_apply_vect.f90 index c55a9f48..af2b926a 100644 --- a/mlprec/impl/smoother/mld_z_as_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_z_as_smoother_apply_vect.f90 @@ -54,12 +54,11 @@ subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& type(psb_z_vect_type),intent(inout), optional :: initu integer(psb_ipk_) :: n_row,n_col, nrow_d, i - complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) - complex(psb_dpk_), allocatable :: vx(:) - type(psb_z_vect_type) :: vtx, vty, vww + complex(psb_dpk_), pointer :: aux(:) + type(psb_z_vect_type) :: tx, ty, ww integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5) character :: trans_, init_ - character(len=20) :: name='z_as_smoother_apply', ch_err + character(len=20) :: name='z_as_smoother_apply_v', ch_err call psb_erractionsave(err_act) @@ -95,55 +94,21 @@ subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& nrow_d = desc_data%get_local_rows() isz = max(n_row,N_COL) - if ((6*isz) <= size(work)) then - ww => work(1:isz) - tx => work(isz+1:2*isz) - ty => work(2*isz+1:3*isz) - aux => work(3*isz+1:) - else if ((4*isz) <= size(work)) then - aux => work(1:) - allocate(ww(isz),tx(isz),ty(isz),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_alloc_request_,name,& - & i_err=(/3*isz,izero,izero,izero,izero/),& - & a_err='complex(psb_dpk_)') - goto 9999 - end if - else if ((3*isz) <= size(work)) then - ww => work(1:isz) - tx => work(isz+1:2*isz) - ty => work(2*isz+1:3*isz) - allocate(aux(4*isz),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_alloc_request_,name,& - & i_err=(/4*isz,izero,izero,izero,izero/),& - & a_err='complex(psb_dpk_)') - goto 9999 - end if + if (4*isz <= size(work)) then + aux => work(:) else - allocate(ww(isz),tx(isz),ty(isz),& - &aux(4*isz),stat=info) + allocate(aux(4*isz),stat=info) if (info /= psb_success_) then call psb_errpush(psb_err_alloc_request_,name,& & i_err=(/4*isz,izero,izero,izero,izero/),& & a_err='complex(psb_dpk_)') goto 9999 end if - endif - if (sweeps == 0) then - + if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then ! - ! K^0 = I - ! zero sweeps of any smoother is just the identity. - ! - call psb_geaxpby(alpha,x,beta,y,desc_data,info) - - else if ((sm%novr == 0).and.(sweeps == 1).and.(.not.sm%sv%is_iterative())) then - ! - ! Shortcut: in this case it's just the same - ! as Block Jacobi. Moreover, if .not.sv%is_iterative, there's no need to pass init + ! Shortcut: in this case there is nothing else to be done. ! call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) @@ -153,388 +118,121 @@ subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& goto 9999 endif - else + else if (sweeps >= 0) then + ! + ! + ! Apply multiple sweeps of an AS solver + ! to compute an approximate solution of a linear system. + ! + ! + call psb_geasb(tx,sm%desc_data,info,mold=x%v,scratch=.true.) + call psb_geasb(ty,sm%desc_data,info,mold=x%v,scratch=.true.) + call psb_geasb(ww,sm%desc_data,info,mold=x%v,scratch=.true.) - call psb_geasb(vtx,sm%desc_data,info,mold=x%v,scratch=.true.) - call psb_geasb(vty,sm%desc_data,info,mold=x%v,scratch=.true.) - call psb_geasb(vww,sm%desc_data,info,mold=x%v,scratch=.true.) + ! + ! Unroll the first iteration and fold it inside SELECT CASE + ! this will save one SPMM when INIT=Z, and will be + ! significant when sweeps=1 (a common case) + ! + call psb_geaxpby(zone,x,zzero,tx,desc_data,info) + if (info == 0) call sm%apply_restr(tx,trans_,aux,info) + if (info == 0) call psb_geaxpby(zone,tx,zzero,ww,sm%desc_data,info) select case (init_) - case('Z') - call vtx%zero() + case('Z') + call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Z') + case('Y') - call psb_geaxpby(zone,y,zzero,vtx,desc_data,info) + call psb_geaxpby(zone,y,zzero,ty,desc_data,info) + if (info == 0) call sm%apply_restr(ty,trans_,aux,info) + if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,& + & work=aux,trans=trans_) + call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,info,init='Y') + case('U') if (.not.present(initu)) then call psb_errpush(psb_err_internal_error_,name,& & a_err='missing initu to smoother_apply') goto 9999 end if - call psb_geaxpby(zone,initu,zzero,vtx,desc_data,info) + call psb_geaxpby(zone,initu,zzero,ty,desc_data,info) + if (info == 0) call sm%apply_restr(ty,trans_,aux,info) + if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,& + & work=aux,trans=trans_) + call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,info,init='Y') + case default call psb_errpush(psb_err_internal_error_,name,& & a_err='wrong init to smoother_apply') goto 9999 end select - - if (sweeps == 1) then - - select case(trans_) - case('N') - ! - ! Get the overlap entries of tx (tx == x) - ! - if (sm%restr == psb_halo_) then - call psb_halo(vtx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - else if (sm%restr /= psb_none_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_restr_') - goto 9999 - end if - - - case('T','C') - ! - ! With transpose, we have to do it here - ! - - select case (sm%prol) - - case(psb_none_) - ! - ! Do nothing - - case(psb_sum_) - ! - ! The transpose of sum is halo - ! - call psb_halo(vtx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - - case(psb_avg_) - ! - ! Tricky one: first we have to scale the overlap entries, - ! which we can do by assignind mode=0, i.e. no communication - ! (hence only scaling), then we do the halo - ! - call psb_ovrl(vtx,sm%desc_data,info,& - & update=psb_avg_,work=aux,mode=izero) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - call psb_halo(vtx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_prol_') - goto 9999 - end select - - - case default - info=psb_err_iarg_invalid_i_ - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - - call sm%sv%apply(zone,vtx,zzero,vty,sm%desc_data,trans_,aux,info,init='Y') - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in sub_aply Jacobi Sweeps = 1') - goto 9999 - endif - - select case(trans_) - case('N') - - select case (sm%prol) - - case(psb_none_) - ! - ! Would work anyway, but since it is supposed to do nothing ... - ! call psb_ovrl(ty,sm%desc_data,info,& - ! & update=sm%prol,work=aux) - - - case(psb_sum_,psb_avg_) - ! - ! Update the overlap of ty - ! - call psb_ovrl(vty,sm%desc_data,info,& - & update=sm%prol,work=aux) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_prol_') - goto 9999 - end select - - case('T','C') - ! - ! With transpose, we have to do it here - ! - if (sm%restr == psb_halo_) then - call psb_ovrl(vty,sm%desc_data,info,& - & update=psb_sum_,work=aux) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - else if (sm%restr /= psb_none_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_restr_') - goto 9999 - end if - - case default - info=psb_err_iarg_invalid_i_ - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - - - - else if (sweeps > 1) then - - ! - ! - ! Apply prec%iprcparm(mld_smoother_sweeps_) sweeps of a block-Jacobi solver - ! to compute an approximate solution of a linear system. + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in sub_aply Jacobi Sweeps = 1') + goto 9999 + endif + + do i=1, sweeps-1 ! + ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the + ! block diagonal part and the remaining part of the local matrix + ! and Y(j) is the approximate solution at sweep j. ! - select case (init_) - case('Z') - call vty%zero() - case('Y') - call psb_geaxpby(zone,y,zzero,vty,sm%desc_data,info) - case('U') - if (.not.present(initu)) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='missing initu to smoother_apply') - goto 9999 - end if - call psb_geaxpby(zone,initu,zzero,vty,sm%desc_data,info) - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='wrong init to smoother_apply') - goto 9999 - end select + if (info == 0) call psb_geaxpby(zone,tx,zzero,ww,sm%desc_data,info) + if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,& + & work=aux,trans=trans_) - do i=1, sweeps - select case(trans_) - case('N') - ! - ! Get the overlap entries of tx (tx == x) - ! - if (sm%restr == psb_halo_) then - call psb_halo(vtx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - else if (sm%restr /= psb_none_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_restr_') - goto 9999 - end if - - - case('T','C') - ! - ! With transpose, we have to do it here - ! - - select case (sm%prol) - - case(psb_none_) - ! - ! Do nothing - - case(psb_sum_) - ! - ! The transpose of sum is halo - ! - call psb_halo(vtx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - - case(psb_avg_) - ! - ! Tricky one: first we have to scale the overlap entries, - ! which we can do by assignind mode=0, i.e. no communication - ! (hence only scaling), then we do the halo - ! - call psb_ovrl(vtx,sm%desc_data,info,& - & update=psb_avg_,work=aux,mode=izero) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - call psb_halo(vtx,sm%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_halo' - goto 9999 - end if - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_prol_') - goto 9999 - end select - - - case default - info=psb_err_iarg_invalid_i_ - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - ! - ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the - ! block diagonal part and the remaining part of the local matrix - ! and Y(j) is the approximate solution at sweep j. - ! - call psb_geaxpby(zone,vtx,zzero,vww,sm%desc_data,info) - call psb_spmm(-zone,sm%nd,vty,zone,vww,sm%desc_data,info,& - & work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call sm%sv%apply(zone,vww,zzero,vty,sm%desc_data,trans_,aux,info,init='Y') - - if (info /= psb_success_) exit - - - select case(trans_) - case('N') - - select case (sm%prol) - - case(psb_none_) - ! - ! Would work anyway, but since it is supposed to do nothing ... - ! call psb_ovrl(ty,sm%desc_data,info,& - ! & update=sm%prol,work=aux) - - - case(psb_sum_,psb_avg_) - ! - ! Update the overlap of ty - ! - call psb_ovrl(vty,sm%desc_data,info,& - & update=sm%prol,work=aux) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_prol_') - goto 9999 - end select + if (info /= psb_success_) exit - case('T','C') - ! - ! With transpose, we have to do it here - ! - if (sm%restr == psb_halo_) then - call psb_ovrl(vty,sm%desc_data,info,& - & update=psb_sum_,work=aux) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_ovrl' - goto 9999 - end if - else if (sm%restr /= psb_none_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Invalid mld_sub_restr_') - goto 9999 - end if - - case default - info=psb_err_iarg_invalid_i_ - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - end do - - if (info /= psb_success_) then - info=psb_err_internal_error_ - call psb_errpush(info,name,& - & a_err='subsolve with Jacobi sweeps > 1') - goto 9999 - end if - - - else + call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Y') + + if (info /= psb_success_) exit + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + + end do - info = psb_err_iarg_neg_ + if (info /= psb_success_) then + info=psb_err_internal_error_ call psb_errpush(info,name,& - & i_err=(/itwo,sweeps,izero,izero,izero/)) - goto 9999 - - + & a_err='subsolve with Jacobi sweeps > 1') + goto 9999 end if - + ! ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) ! - call psb_geaxpby(alpha,vty,beta,y,desc_data,info) - - end if + call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + + + else + + info = psb_err_iarg_neg_ + call psb_errpush(info,name,& + & i_err=(/itwo,sweeps,izero,izero,izero/)) + goto 9999 + endif - if ((6*isz) <= size(work)) then - else if ((4*isz) <= size(work)) then - deallocate(ww,tx,ty) - else if ((3*isz) <= size(work)) then - deallocate(aux) - else - deallocate(ww,aux,tx,ty) + + if (.not.(4*isz <= size(work))) then + deallocate(aux,stat=info) endif - call vww%free(info) - call vtx%free(info) - call vty%free(info) + if (info ==0) call ww%free(info) + if (info ==0) call tx%free(info) + if (info ==0) call ty%free(info) + if (info /= 0) then + info = psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + end if call psb_erractionrestore(err_act) return - + 9999 call psb_error_handler(err_act) - + return end subroutine mld_z_as_smoother_apply_vect diff --git a/mlprec/impl/smoother/mld_z_as_smoother_prol_a.f90 b/mlprec/impl/smoother/mld_z_as_smoother_prol_a.f90 new file mode 100644 index 00000000..c0118e98 --- /dev/null +++ b/mlprec/impl/smoother/mld_z_as_smoother_prol_a.f90 @@ -0,0 +1,150 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ 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_as_smoother_prol_a(sm,x,trans,work,info,data) + use psb_base_mod + use mld_z_as_smoother, mld_protect_nam => mld_z_as_smoother_prol_a + implicit none + class(mld_z_as_smoother_type), intent(inout) :: sm + complex(psb_dpk_), intent(inout) :: x(:) + character(len=1),intent(in) :: trans + complex(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: data + !Local + integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5), data_ + character :: trans_ + character(len=20) :: name='z_as_smther_prol_a', ch_err + + call psb_erractionsave(err_act) + + info = psb_success_ + ictxt = sm%desc_data%get_context() + call psb_info(ictxt,me,np) + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T') + case('C') + case default + info = psb_err_iarg_invalid_i_ + call psb_errpush(info,name) + goto 9999 + end select + + if (present(data)) then + data_ = data + else + data_ = psb_comm_ext_ + end if + + if (.not.allocated(sm%sv)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + + + select case(trans_) + case('N') + + select case (sm%prol) + + case(psb_none_) + ! + ! Would work anyway, but since it is supposed to do nothing ... + ! call psb_ovrl(x,sm%desc_data,info,& + ! & update=sm%prol,work=work) + + + case(psb_sum_,psb_avg_) + ! + ! Update the overlap of x + ! + call psb_ovrl(x,sm%desc_data,info,& + & update=sm%prol,work=work) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_ovrl' + goto 9999 + end if + + case default + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Invalid mld_sub_prol_') + goto 9999 + end select + + case('T','C') + ! + ! With transpose, we have to do it here + ! + if (sm%restr == psb_halo_) then + call psb_ovrl(x,sm%desc_data,info,& + & update=psb_sum_,work=work) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_ovrl' + goto 9999 + end if + else if (sm%restr /= psb_none_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Invalid mld_sub_restr_') + goto 9999 + end if + + case default + info=psb_err_iarg_invalid_i_ + int_err(1)=6 + ch_err(2:2)=trans + goto 9999 + end select + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_z_as_smoother_prol_a + + diff --git a/mlprec/impl/smoother/mld_z_as_smoother_prol_v.f90 b/mlprec/impl/smoother/mld_z_as_smoother_prol_v.f90 new file mode 100644 index 00000000..c8bdec78 --- /dev/null +++ b/mlprec/impl/smoother/mld_z_as_smoother_prol_v.f90 @@ -0,0 +1,150 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ 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_as_smoother_prol_v(sm,x,trans,work,info,data) + use psb_base_mod + use mld_z_as_smoother, mld_protect_nam => mld_z_as_smoother_prol_v + implicit none + class(mld_z_as_smoother_type), intent(inout) :: sm + type(psb_z_vect_type),intent(inout) :: x + character(len=1),intent(in) :: trans + complex(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: data + !Local + integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5), data_ + character :: trans_ + character(len=20) :: name='z_as_smther_prol_v', ch_err + + call psb_erractionsave(err_act) + + info = psb_success_ + ictxt = sm%desc_data%get_context() + call psb_info(ictxt,me,np) + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T') + case('C') + case default + info = psb_err_iarg_invalid_i_ + call psb_errpush(info,name) + goto 9999 + end select + + if (present(data)) then + data_ = data + else + data_ = psb_comm_ext_ + end if + + if (.not.allocated(sm%sv)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + + + select case(trans_) + case('N') + + select case (sm%prol) + + case(psb_none_) + ! + ! Would work anyway, but since it is supposed to do nothing ... + ! call psb_ovrl(x,sm%desc_data,info,& + ! & update=sm%prol,work=work) + + + case(psb_sum_,psb_avg_) + ! + ! Update the overlap of x + ! + call psb_ovrl(x,sm%desc_data,info,& + & update=sm%prol,work=work) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_ovrl' + goto 9999 + end if + + case default + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Invalid mld_sub_prol_') + goto 9999 + end select + + case('T','C') + ! + ! With transpose, we have to do it here + ! + if (sm%restr == psb_halo_) then + call psb_ovrl(x,sm%desc_data,info,& + & update=psb_sum_,work=work) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_ovrl' + goto 9999 + end if + else if (sm%restr /= psb_none_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Invalid mld_sub_restr_') + goto 9999 + end if + + case default + info=psb_err_iarg_invalid_i_ + int_err(1)=6 + ch_err(2:2)=trans + goto 9999 + end select + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_z_as_smoother_prol_v + + diff --git a/mlprec/impl/smoother/mld_z_as_smoother_restr_a.f90 b/mlprec/impl/smoother/mld_z_as_smoother_restr_a.f90 new file mode 100644 index 00000000..a7259144 --- /dev/null +++ b/mlprec/impl/smoother/mld_z_as_smoother_restr_a.f90 @@ -0,0 +1,169 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ 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_as_smoother_restr_a(sm,x,trans,work,info,data) + use psb_base_mod + use mld_z_as_smoother, mld_protect_nam => mld_z_as_smoother_restr_a + implicit none + class(mld_z_as_smoother_type), intent(inout) :: sm + complex(psb_dpk_), intent(inout) :: x(:) + character(len=1),intent(in) :: trans + complex(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: data + !Local + integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5), data_ + character :: trans_ + character(len=20) :: name='z_as_smther_restr_a', ch_err + + call psb_erractionsave(err_act) + + info = psb_success_ + ictxt = sm%desc_data%get_context() + call psb_info(ictxt,me,np) + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T') + case('C') + case default + info = psb_err_iarg_invalid_i_ + call psb_errpush(info,name) + goto 9999 + end select + + if (present(data)) then + data_ = data + else + data_ = psb_comm_ext_ + end if + + if (.not.allocated(sm%sv)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + + select case(trans_) + case('N') + ! + ! Get the overlap entries x + ! + if (sm%restr == psb_halo_) then + call psb_halo(x,sm%desc_data,info,work=work,data=data_) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_halo' + goto 9999 + end if + else if (sm%restr /= psb_none_) then + call psb_errpush(psb_err_internal_error_,name,& + &a_err='Invalid mld_sub_restr_') + goto 9999 + end if + + + case('T','C') + ! + ! With transpose, we have to do it here + ! + + select case (sm%prol) + + case(psb_none_) + ! + ! Do nothing + + case(psb_sum_) + ! + ! The transpose of sum is halo + ! + call psb_halo(x,sm%desc_data,info,work=work,data=data_) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_halo' + goto 9999 + end if + + case(psb_avg_) + ! + ! Tricky one: first we have to scale the overlap entries, + ! which we can do by assignind mode=0, i.e. no communication + ! (hence only scaling), then we do the halo + ! + call psb_ovrl(x,sm%desc_data,info,& + & update=psb_avg_,work=work,mode=izero) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_ovrl' + goto 9999 + end if + call psb_halo(x,sm%desc_data,info,work=work,data=data_) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_halo' + goto 9999 + end if + + case default + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Invalid mld_sub_prol_') + goto 9999 + end select + + + case default + info=psb_err_iarg_invalid_i_ + int_err(1)=6 + ch_err(2:2)=trans + goto 9999 + end select + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_z_as_smoother_restr_a + + diff --git a/mlprec/impl/smoother/mld_z_as_smoother_restr_v.f90 b/mlprec/impl/smoother/mld_z_as_smoother_restr_v.f90 new file mode 100644 index 00000000..0424830f --- /dev/null +++ b/mlprec/impl/smoother/mld_z_as_smoother_restr_v.f90 @@ -0,0 +1,169 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ 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_as_smoother_restr_v(sm,x,trans,work,info,data) + use psb_base_mod + use mld_z_as_smoother, mld_protect_nam => mld_z_as_smoother_restr_v + implicit none + class(mld_z_as_smoother_type), intent(inout) :: sm + type(psb_z_vect_type),intent(inout) :: x + character(len=1),intent(in) :: trans + complex(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: data + !Local + integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5), data_ + character :: trans_ + character(len=20) :: name='z_as_smther_restr_v', ch_err + + call psb_erractionsave(err_act) + + info = psb_success_ + ictxt = sm%desc_data%get_context() + call psb_info(ictxt,me,np) + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T') + case('C') + case default + info = psb_err_iarg_invalid_i_ + call psb_errpush(info,name) + goto 9999 + end select + + if (present(data)) then + data_ = data + else + data_ = psb_comm_ext_ + end if + + if (.not.allocated(sm%sv)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + + select case(trans_) + case('N') + ! + ! Get the overlap entries x + ! + if (sm%restr == psb_halo_) then + call psb_halo(x,sm%desc_data,info,work=work,data=data_) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_halo' + goto 9999 + end if + else if (sm%restr /= psb_none_) then + call psb_errpush(psb_err_internal_error_,name,& + &a_err='Invalid mld_sub_restr_') + goto 9999 + end if + + + case('T','C') + ! + ! With transpose, we have to do it here + ! + + select case (sm%prol) + + case(psb_none_) + ! + ! Do nothing + + case(psb_sum_) + ! + ! The transpose of sum is halo + ! + call psb_halo(x,sm%desc_data,info,work=work,data=data_) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_halo' + goto 9999 + end if + + case(psb_avg_) + ! + ! Tricky one: first we have to scale the overlap entries, + ! which we can do by assignind mode=0, i.e. no communication + ! (hence only scaling), then we do the halo + ! + call psb_ovrl(x,sm%desc_data,info,& + & update=psb_avg_,work=work,mode=izero) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_ovrl' + goto 9999 + end if + call psb_halo(x,sm%desc_data,info,work=work,data=data_) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_halo' + goto 9999 + end if + + case default + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Invalid mld_sub_prol_') + goto 9999 + end select + + + case default + info=psb_err_iarg_invalid_i_ + int_err(1)=6 + ch_err(2:2)=trans + goto 9999 + end select + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_z_as_smoother_restr_v + + diff --git a/mlprec/mld_c_as_smoother.f90 b/mlprec/mld_c_as_smoother.f90 index 9c8bcf5a..9b968a00 100644 --- a/mlprec/mld_c_as_smoother.f90 +++ b/mlprec/mld_c_as_smoother.f90 @@ -62,6 +62,12 @@ module mld_c_as_smoother procedure, pass(sm) :: clone => mld_c_as_smoother_clone procedure, pass(sm) :: apply_v => mld_c_as_smoother_apply_vect procedure, pass(sm) :: apply_a => mld_c_as_smoother_apply + procedure, pass(sm) :: restr_a => mld_c_as_smoother_restr_a + procedure, pass(sm) :: prol_a => mld_c_as_smoother_prol_a + procedure, pass(sm) :: restr_v => mld_c_as_smoother_restr_v + procedure, pass(sm) :: prol_v => mld_c_as_smoother_prol_v + generic, public :: apply_restr => restr_v, restr_a + generic, public :: apply_prol => prol_v, prol_a procedure, pass(sm) :: free => mld_c_as_smoother_free procedure, pass(sm) :: seti => mld_c_as_smoother_seti procedure, pass(sm) :: setc => mld_c_as_smoother_setc @@ -95,6 +101,67 @@ module mld_c_as_smoother end subroutine mld_c_as_smoother_check end interface + interface + subroutine mld_c_as_smoother_restr_v(sm,x,trans,work,info,data) + import :: psb_cspmat_type, psb_c_vect_type, psb_c_base_vect_type, & + & psb_spk_, mld_c_as_smoother_type, psb_long_int_k_, & + & psb_desc_type, psb_ipk_ + implicit none + class(mld_c_as_smoother_type), intent(inout) :: sm + type(psb_c_vect_type),intent(inout) :: x + character(len=1),intent(in) :: trans + complex(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: data + end subroutine mld_c_as_smoother_restr_v + end interface + + interface + subroutine mld_c_as_smoother_restr_a(sm,x,trans,work,info,data) + import :: psb_cspmat_type, psb_c_vect_type, psb_c_base_vect_type, & + & psb_spk_, mld_c_as_smoother_type, psb_long_int_k_, & + & psb_desc_type, psb_ipk_ + implicit none + class(mld_c_as_smoother_type), intent(inout) :: sm + complex(psb_spk_), intent(inout) :: x(:) + character(len=1),intent(in) :: trans + complex(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: data + end subroutine mld_c_as_smoother_restr_a + end interface + + interface + subroutine mld_c_as_smoother_prol_v(sm,x,trans,work,info,data) + import :: psb_cspmat_type, psb_c_vect_type, psb_c_base_vect_type, & + & psb_spk_, mld_c_as_smoother_type, psb_long_int_k_, & + & psb_desc_type, psb_ipk_ + implicit none + class(mld_c_as_smoother_type), intent(inout) :: sm + type(psb_c_vect_type),intent(inout) :: x + character(len=1),intent(in) :: trans + complex(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: data + end subroutine mld_c_as_smoother_prol_v + end interface + + interface + subroutine mld_c_as_smoother_prol_a(sm,x,trans,work,info,data) + import :: psb_cspmat_type, psb_c_vect_type, psb_c_base_vect_type, & + & psb_spk_, mld_c_as_smoother_type, psb_long_int_k_, & + & psb_desc_type, psb_ipk_ + implicit none + class(mld_c_as_smoother_type), intent(inout) :: sm + complex(psb_spk_), intent(inout) :: x(:) + character(len=1),intent(in) :: trans + complex(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: data + end subroutine mld_c_as_smoother_prol_a + end interface + + interface subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& & trans,sweeps,work,info,init,initu) diff --git a/mlprec/mld_d_as_smoother.f90 b/mlprec/mld_d_as_smoother.f90 index f67d3576..b6282334 100644 --- a/mlprec/mld_d_as_smoother.f90 +++ b/mlprec/mld_d_as_smoother.f90 @@ -62,6 +62,12 @@ module mld_d_as_smoother procedure, pass(sm) :: clone => mld_d_as_smoother_clone procedure, pass(sm) :: apply_v => mld_d_as_smoother_apply_vect procedure, pass(sm) :: apply_a => mld_d_as_smoother_apply + procedure, pass(sm) :: restr_a => mld_d_as_smoother_restr_a + procedure, pass(sm) :: prol_a => mld_d_as_smoother_prol_a + procedure, pass(sm) :: restr_v => mld_d_as_smoother_restr_v + procedure, pass(sm) :: prol_v => mld_d_as_smoother_prol_v + generic, public :: apply_restr => restr_v, restr_a + generic, public :: apply_prol => prol_v, prol_a procedure, pass(sm) :: free => mld_d_as_smoother_free procedure, pass(sm) :: seti => mld_d_as_smoother_seti procedure, pass(sm) :: setc => mld_d_as_smoother_setc @@ -95,6 +101,67 @@ module mld_d_as_smoother end subroutine mld_d_as_smoother_check end interface + interface + subroutine mld_d_as_smoother_restr_v(sm,x,trans,work,info,data) + import :: psb_dspmat_type, psb_d_vect_type, psb_d_base_vect_type, & + & psb_dpk_, mld_d_as_smoother_type, psb_long_int_k_, & + & psb_desc_type, psb_ipk_ + implicit none + class(mld_d_as_smoother_type), intent(inout) :: sm + type(psb_d_vect_type),intent(inout) :: x + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: data + end subroutine mld_d_as_smoother_restr_v + end interface + + interface + subroutine mld_d_as_smoother_restr_a(sm,x,trans,work,info,data) + import :: psb_dspmat_type, psb_d_vect_type, psb_d_base_vect_type, & + & psb_dpk_, mld_d_as_smoother_type, psb_long_int_k_, & + & psb_desc_type, psb_ipk_ + implicit none + class(mld_d_as_smoother_type), intent(inout) :: sm + real(psb_dpk_), intent(inout) :: x(:) + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: data + end subroutine mld_d_as_smoother_restr_a + end interface + + interface + subroutine mld_d_as_smoother_prol_v(sm,x,trans,work,info,data) + import :: psb_dspmat_type, psb_d_vect_type, psb_d_base_vect_type, & + & psb_dpk_, mld_d_as_smoother_type, psb_long_int_k_, & + & psb_desc_type, psb_ipk_ + implicit none + class(mld_d_as_smoother_type), intent(inout) :: sm + type(psb_d_vect_type),intent(inout) :: x + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: data + end subroutine mld_d_as_smoother_prol_v + end interface + + interface + subroutine mld_d_as_smoother_prol_a(sm,x,trans,work,info,data) + import :: psb_dspmat_type, psb_d_vect_type, psb_d_base_vect_type, & + & psb_dpk_, mld_d_as_smoother_type, psb_long_int_k_, & + & psb_desc_type, psb_ipk_ + implicit none + class(mld_d_as_smoother_type), intent(inout) :: sm + real(psb_dpk_), intent(inout) :: x(:) + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: data + end subroutine mld_d_as_smoother_prol_a + end interface + + interface subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& & trans,sweeps,work,info,init,initu) diff --git a/mlprec/mld_s_as_smoother.f90 b/mlprec/mld_s_as_smoother.f90 index 82d04e38..17b010e8 100644 --- a/mlprec/mld_s_as_smoother.f90 +++ b/mlprec/mld_s_as_smoother.f90 @@ -62,6 +62,12 @@ module mld_s_as_smoother procedure, pass(sm) :: clone => mld_s_as_smoother_clone procedure, pass(sm) :: apply_v => mld_s_as_smoother_apply_vect procedure, pass(sm) :: apply_a => mld_s_as_smoother_apply + procedure, pass(sm) :: restr_a => mld_s_as_smoother_restr_a + procedure, pass(sm) :: prol_a => mld_s_as_smoother_prol_a + procedure, pass(sm) :: restr_v => mld_s_as_smoother_restr_v + procedure, pass(sm) :: prol_v => mld_s_as_smoother_prol_v + generic, public :: apply_restr => restr_v, restr_a + generic, public :: apply_prol => prol_v, prol_a procedure, pass(sm) :: free => mld_s_as_smoother_free procedure, pass(sm) :: seti => mld_s_as_smoother_seti procedure, pass(sm) :: setc => mld_s_as_smoother_setc @@ -95,6 +101,67 @@ module mld_s_as_smoother end subroutine mld_s_as_smoother_check end interface + interface + subroutine mld_s_as_smoother_restr_v(sm,x,trans,work,info,data) + import :: psb_sspmat_type, psb_s_vect_type, psb_s_base_vect_type, & + & psb_spk_, mld_s_as_smoother_type, psb_long_int_k_, & + & psb_desc_type, psb_ipk_ + implicit none + class(mld_s_as_smoother_type), intent(inout) :: sm + type(psb_s_vect_type),intent(inout) :: x + character(len=1),intent(in) :: trans + real(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: data + end subroutine mld_s_as_smoother_restr_v + end interface + + interface + subroutine mld_s_as_smoother_restr_a(sm,x,trans,work,info,data) + import :: psb_sspmat_type, psb_s_vect_type, psb_s_base_vect_type, & + & psb_spk_, mld_s_as_smoother_type, psb_long_int_k_, & + & psb_desc_type, psb_ipk_ + implicit none + class(mld_s_as_smoother_type), intent(inout) :: sm + real(psb_spk_), intent(inout) :: x(:) + character(len=1),intent(in) :: trans + real(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: data + end subroutine mld_s_as_smoother_restr_a + end interface + + interface + subroutine mld_s_as_smoother_prol_v(sm,x,trans,work,info,data) + import :: psb_sspmat_type, psb_s_vect_type, psb_s_base_vect_type, & + & psb_spk_, mld_s_as_smoother_type, psb_long_int_k_, & + & psb_desc_type, psb_ipk_ + implicit none + class(mld_s_as_smoother_type), intent(inout) :: sm + type(psb_s_vect_type),intent(inout) :: x + character(len=1),intent(in) :: trans + real(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: data + end subroutine mld_s_as_smoother_prol_v + end interface + + interface + subroutine mld_s_as_smoother_prol_a(sm,x,trans,work,info,data) + import :: psb_sspmat_type, psb_s_vect_type, psb_s_base_vect_type, & + & psb_spk_, mld_s_as_smoother_type, psb_long_int_k_, & + & psb_desc_type, psb_ipk_ + implicit none + class(mld_s_as_smoother_type), intent(inout) :: sm + real(psb_spk_), intent(inout) :: x(:) + character(len=1),intent(in) :: trans + real(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: data + end subroutine mld_s_as_smoother_prol_a + end interface + + interface subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& & trans,sweeps,work,info,init,initu) diff --git a/mlprec/mld_z_as_smoother.f90 b/mlprec/mld_z_as_smoother.f90 index 0aa57cdb..cfe6ee78 100644 --- a/mlprec/mld_z_as_smoother.f90 +++ b/mlprec/mld_z_as_smoother.f90 @@ -62,6 +62,12 @@ module mld_z_as_smoother procedure, pass(sm) :: clone => mld_z_as_smoother_clone procedure, pass(sm) :: apply_v => mld_z_as_smoother_apply_vect procedure, pass(sm) :: apply_a => mld_z_as_smoother_apply + procedure, pass(sm) :: restr_a => mld_z_as_smoother_restr_a + procedure, pass(sm) :: prol_a => mld_z_as_smoother_prol_a + procedure, pass(sm) :: restr_v => mld_z_as_smoother_restr_v + procedure, pass(sm) :: prol_v => mld_z_as_smoother_prol_v + generic, public :: apply_restr => restr_v, restr_a + generic, public :: apply_prol => prol_v, prol_a procedure, pass(sm) :: free => mld_z_as_smoother_free procedure, pass(sm) :: seti => mld_z_as_smoother_seti procedure, pass(sm) :: setc => mld_z_as_smoother_setc @@ -95,6 +101,67 @@ module mld_z_as_smoother end subroutine mld_z_as_smoother_check end interface + interface + subroutine mld_z_as_smoother_restr_v(sm,x,trans,work,info,data) + import :: psb_zspmat_type, psb_z_vect_type, psb_z_base_vect_type, & + & psb_dpk_, mld_z_as_smoother_type, psb_long_int_k_, & + & psb_desc_type, psb_ipk_ + implicit none + class(mld_z_as_smoother_type), intent(inout) :: sm + type(psb_z_vect_type),intent(inout) :: x + character(len=1),intent(in) :: trans + complex(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: data + end subroutine mld_z_as_smoother_restr_v + end interface + + interface + subroutine mld_z_as_smoother_restr_a(sm,x,trans,work,info,data) + import :: psb_zspmat_type, psb_z_vect_type, psb_z_base_vect_type, & + & psb_dpk_, mld_z_as_smoother_type, psb_long_int_k_, & + & psb_desc_type, psb_ipk_ + implicit none + class(mld_z_as_smoother_type), intent(inout) :: sm + complex(psb_dpk_), intent(inout) :: x(:) + character(len=1),intent(in) :: trans + complex(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: data + end subroutine mld_z_as_smoother_restr_a + end interface + + interface + subroutine mld_z_as_smoother_prol_v(sm,x,trans,work,info,data) + import :: psb_zspmat_type, psb_z_vect_type, psb_z_base_vect_type, & + & psb_dpk_, mld_z_as_smoother_type, psb_long_int_k_, & + & psb_desc_type, psb_ipk_ + implicit none + class(mld_z_as_smoother_type), intent(inout) :: sm + type(psb_z_vect_type),intent(inout) :: x + character(len=1),intent(in) :: trans + complex(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: data + end subroutine mld_z_as_smoother_prol_v + end interface + + interface + subroutine mld_z_as_smoother_prol_a(sm,x,trans,work,info,data) + import :: psb_zspmat_type, psb_z_vect_type, psb_z_base_vect_type, & + & psb_dpk_, mld_z_as_smoother_type, psb_long_int_k_, & + & psb_desc_type, psb_ipk_ + implicit none + class(mld_z_as_smoother_type), intent(inout) :: sm + complex(psb_dpk_), intent(inout) :: x(:) + character(len=1),intent(in) :: trans + complex(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), optional, intent(in) :: data + end subroutine mld_z_as_smoother_prol_a + end interface + + interface subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& & trans,sweeps,work,info,init,initu)