From c9e17190988182aacdc3b782fe6a950648714072 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Tue, 9 Jun 2020 10:23:37 +0200 Subject: [PATCH] First round of changes for Richardson implementation. Zero sweeps is the identity. --- .../impl/smoother/mld_c_as_smoother_bld.f90 | 2 +- .../impl/smoother/mld_c_as_smoother_clone.f90 | 1 + .../impl/smoother/mld_c_as_smoother_free.f90 | 1 + .../smoother/mld_c_jac_smoother_apply.f90 | 179 ++++---------- .../mld_c_jac_smoother_apply_vect.f90 | 234 +++++------------- .../impl/smoother/mld_c_jac_smoother_bld.f90 | 5 +- .../smoother/mld_c_jac_smoother_clone.f90 | 1 + .../impl/smoother/mld_d_as_smoother_bld.f90 | 2 +- .../impl/smoother/mld_d_as_smoother_clone.f90 | 1 + .../impl/smoother/mld_d_as_smoother_free.f90 | 1 + .../smoother/mld_d_jac_smoother_apply.f90 | 179 ++++---------- .../mld_d_jac_smoother_apply_vect.f90 | 234 +++++------------- .../impl/smoother/mld_d_jac_smoother_bld.f90 | 5 +- .../smoother/mld_d_jac_smoother_clone.f90 | 1 + .../impl/smoother/mld_s_as_smoother_bld.f90 | 2 +- .../impl/smoother/mld_s_as_smoother_clone.f90 | 1 + .../impl/smoother/mld_s_as_smoother_free.f90 | 1 + .../smoother/mld_s_jac_smoother_apply.f90 | 179 ++++---------- .../mld_s_jac_smoother_apply_vect.f90 | 234 +++++------------- .../impl/smoother/mld_s_jac_smoother_bld.f90 | 5 +- .../smoother/mld_s_jac_smoother_clone.f90 | 1 + .../impl/smoother/mld_z_as_smoother_bld.f90 | 2 +- .../impl/smoother/mld_z_as_smoother_clone.f90 | 1 + .../impl/smoother/mld_z_as_smoother_free.f90 | 1 + .../smoother/mld_z_jac_smoother_apply.f90 | 179 ++++---------- .../mld_z_jac_smoother_apply_vect.f90 | 234 +++++------------- .../impl/smoother/mld_z_jac_smoother_bld.f90 | 5 +- .../smoother/mld_z_jac_smoother_clone.f90 | 1 + mlprec/mld_c_as_smoother.f90 | 1 + mlprec/mld_d_as_smoother.f90 | 1 + mlprec/mld_s_as_smoother.f90 | 1 + mlprec/mld_z_as_smoother.f90 | 1 + 32 files changed, 500 insertions(+), 1196 deletions(-) diff --git a/mlprec/impl/smoother/mld_c_as_smoother_bld.f90 b/mlprec/impl/smoother/mld_c_as_smoother_bld.f90 index 19dff29d..e5d17961 100644 --- a/mlprec/impl/smoother/mld_c_as_smoother_bld.f90 +++ b/mlprec/impl/smoother/mld_c_as_smoother_bld.f90 @@ -66,7 +66,7 @@ subroutine mld_c_as_smoother_bld(a,desc_a,sm,info,amold,vmold,imold) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),' start' - + sm%pa => a novr = sm%novr if (novr < 0) then info=psb_err_invalid_ovr_num_ diff --git a/mlprec/impl/smoother/mld_c_as_smoother_clone.f90 b/mlprec/impl/smoother/mld_c_as_smoother_clone.f90 index badcf5ca..572ae031 100644 --- a/mlprec/impl/smoother/mld_c_as_smoother_clone.f90 +++ b/mlprec/impl/smoother/mld_c_as_smoother_clone.f90 @@ -77,6 +77,7 @@ subroutine mld_c_as_smoother_clone(sm,smout,info) allocate(smout%sv,mold=sm%sv,stat=info) if (info == psb_success_) call sm%sv%clone(smo%sv,info) end if + smo%pa => sm%pa class default info = psb_err_internal_error_ diff --git a/mlprec/impl/smoother/mld_c_as_smoother_free.f90 b/mlprec/impl/smoother/mld_c_as_smoother_free.f90 index 3f299a89..51a6b1a5 100644 --- a/mlprec/impl/smoother/mld_c_as_smoother_free.f90 +++ b/mlprec/impl/smoother/mld_c_as_smoother_free.f90 @@ -61,6 +61,7 @@ subroutine mld_c_as_smoother_free(sm,info) end if end if call sm%nd%free() + sm%pa => null() call psb_erractionrestore(err_act) return diff --git a/mlprec/impl/smoother/mld_c_jac_smoother_apply.f90 b/mlprec/impl/smoother/mld_c_jac_smoother_apply.f90 index e76c6707..c7ae272a 100644 --- a/mlprec/impl/smoother/mld_c_jac_smoother_apply.f90 +++ b/mlprec/impl/smoother/mld_c_jac_smoother_apply.f90 @@ -65,7 +65,7 @@ subroutine mld_c_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,& ictxt = desc_data%get_context() call psb_info(ictxt,me,np) - + if (present(init)) then init_ = psb_toupper(init) else @@ -102,146 +102,63 @@ subroutine mld_c_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,& goto 9999 end if endif - - if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then - ! if .not.sv%is_iterative, there's no need to pass init - call sm%sv%apply(alpha,x,beta,y,desc_data,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 - - else if (sweeps >= 0) then - if (associated(sm%pa)) then - ! - ! This means we are dealing with a pure Jacobi smoother/solver. - ! - call psb_geasb(tx,desc_data,info) - call psb_geasb(ty,desc_data,info) - select case (init_) - case('Z') - - call sm%sv%apply(cone,x,czero,ty,desc_data,trans_,aux,info,init='Z') - - case('Y') - call psb_geaxpby(cone,x,czero,tx,desc_data,info) - call psb_geaxpby(cone,y,czero,ty,desc_data,info) - call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(cone,tx,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,x,czero,tx,desc_data,info) - call psb_geaxpby(cone,initu,czero,ty,desc_data,info) - call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,info,init='Y') - - case default + + if (sweeps > 0) then + call psb_geasb(tx,desc_data,info) + call psb_geasb(ty,desc_data,info) + select case (init_) + case('Z') + + call sm%sv%apply(cone,x,czero,ty,desc_data,trans_,aux,info,init='Z') + + case('Y') + call psb_geaxpby(cone,x,czero,tx,desc_data,info) + call psb_geaxpby(cone,y,czero,ty,desc_data,info) + call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(cone,tx,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='wrong init to smoother_apply') + & a_err='missing initu to smoother_apply') goto 9999 - end select - - do i=1, sweeps-1 - ! - ! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)), - ! where is the diagonal and A the matrix. - ! - call psb_geaxpby(cone,x,czero,tx,desc_data,info) - call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call sm%sv%apply(cone,tx,cone,ty,desc_data,trans_,aux,info,init='Y') - - if (info /= psb_success_) exit - end do - - - if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) - - 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 - ! - ! - ! Apply multiple sweeps of a block-Jacobi solver - ! to compute an approximate solution of a linear system. - ! - ! - call psb_geasb(tx,desc_data,info) - call psb_geasb(ty,desc_data,info) + call psb_geaxpby(cone,x,czero,tx,desc_data,info) + call psb_geaxpby(cone,initu,czero,ty,desc_data,info) + call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(cone,tx,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 + do i=1, sweeps-1 ! - ! Unroll the first iteration and fold it inside SELECT CASE - ! this will save one AXPBY and one SPMM when INIT=Z, and will be - ! significant when sweeps=1 (a common case) + ! Compute Y(j+1) = Y(j)+ M^(-1)*(X-A*Y(j)), ! - select case (init_) - case('Z') - - call sm%sv%apply(cone,x,czero,ty,desc_data,trans_,aux,info,init='Z') - - case('Y') - call psb_geaxpby(cone,x,czero,tx,desc_data,info) - call psb_geaxpby(cone,y,czero,ty,desc_data,info) - call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(cone,tx,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,x,czero,tx,desc_data,info) - call psb_geaxpby(cone,initu,czero,ty,desc_data,info) - call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(cone,tx,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 + call psb_geaxpby(cone,x,czero,tx,desc_data,info) + call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_) - 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. - ! - call psb_geaxpby(cone,x,czero,tx,desc_data,info) - call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_) + if (info /= psb_success_) exit - if (info /= psb_success_) exit + call sm%sv%apply(cone,tx,cone,ty,desc_data,trans_,aux,info,init='Y') - call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,info,init='Y') + if (info /= psb_success_) exit + end do - if (info /= psb_success_) exit - end do - if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) - 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 + 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 - + + deallocate(tx,ty,stat=info) if (info /= psb_success_) then info=psb_err_internal_error_ @@ -250,6 +167,12 @@ subroutine mld_c_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,& goto 9999 end if + else if (sweeps == 0) then + ! + ! 0 sweeps of smoother is the identity operator + ! + call psb_geaxpby(alpha,x,beta,y,desc_data,info) + else info = psb_err_iarg_neg_ diff --git a/mlprec/impl/smoother/mld_c_jac_smoother_apply_vect.f90 b/mlprec/impl/smoother/mld_c_jac_smoother_apply_vect.f90 index 2b55f98d..2b9ae54a 100644 --- a/mlprec/impl/smoother/mld_c_jac_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_c_jac_smoother_apply_vect.f90 @@ -49,7 +49,7 @@ subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& type(psb_c_vect_type),intent(inout) :: y complex(psb_spk_),intent(in) :: alpha,beta character(len=1),intent(in) :: trans - integer(psb_ipk_), intent(in) :: sweeps + integer(psb_ipk_), intent(in) :: sweeps complex(psb_spk_),target, intent(inout) :: work(:) type(psb_c_vect_type),intent(inout) :: wv(:) integer(psb_ipk_), intent(out) :: info @@ -108,190 +108,90 @@ subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& end if endif - if(sm%checkres) then - call psb_geall(r,desc_data,info) - call psb_geasb(r,desc_data,info) - resdenum = psb_genrm2(x,desc_data,info) - end if + if (sweeps > 0) then - if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then - ! if .not.sv%is_iterative, there's no need to pass init - call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info) + if(sm%checkres) then + call psb_geall(r,desc_data,info) + call psb_geasb(r,desc_data,info) + resdenum = psb_genrm2(x,desc_data,info) + end if - 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 - - else if (sweeps >= 0) then - select type (smsv => sm%sv) - class is (mld_c_diag_solver_type) - ! - ! This means we are dealing with a pure Jacobi smoother/solver. - ! - associate(tx => wv(1), ty => wv(2)) - select case (init_) - case('Z') - - call sm%sv%apply(cone,x,czero,ty,desc_data,trans_,aux,wv(3:),info,init='Z') - - case('Y') - call psb_geaxpby(cone,x,czero,tx,desc_data,info) - call psb_geaxpby(cone,y,czero,ty,desc_data,info) - call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,wv(3:),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,x,czero,tx,desc_data,info) - call psb_geaxpby(cone,initu,czero,ty,desc_data,info) - call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') + associate(tx => wv(1), ty => wv(2)) + select case (init_) + case('Z') + + call sm%sv%apply(cone,x,czero,ty,desc_data,trans_,aux,wv(3:),info,init='Z') + + case('Y') + call psb_geaxpby(cone,x,czero,tx,desc_data,info) + call psb_geaxpby(cone,y,czero,ty,desc_data,info) + call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') - case default + case('U') + if (.not.present(initu)) then call psb_errpush(psb_err_internal_error_,name,& - & a_err='wrong init to smoother_apply') + & a_err='missing initu to smoother_apply') goto 9999 - end select - - do i=1, sweeps-1 - ! - ! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)), - ! where is the diagonal and A the matrix. - ! - call psb_geaxpby(cone,x,czero,tx,desc_data,info) - call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call sm%sv%apply(cone,tx,cone,ty,desc_data,trans_,aux,wv(3:),info,init='Y') - - if (info /= psb_success_) exit - - if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then - call psb_geaxpby(cone,x,czero,r,r,desc_data,info) - call psb_spmm(-cone,sm%pa,ty,cone,r,desc_data,info) - res = psb_genrm2(r,desc_data,info) - if( sm%printres ) then - call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol) - end if - if ( res < sm%tol*resdenum ) then - if( (sm%printres).and.(mod(sm%printiter,sm%checkiter)/=0) ) & - & call log_conv("BJAC",me,i,1,res,resdenum,sm%tol) - exit - end if - end if + end if + call psb_geaxpby(cone,x,czero,tx,desc_data,info) + call psb_geaxpby(cone,initu,czero,ty,desc_data,info) + call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') + + case default + call psb_errpush(psb_err_internal_error_,name,& + & a_err='wrong init to smoother_apply') + goto 9999 + end select - end do + do i=1, sweeps-1 + ! + ! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)), + ! where is the diagonal and A the matrix. + ! + call psb_geaxpby(cone,x,czero,tx,desc_data,info) + call psb_spmm(-cone,sm%pa,ty,cone,tx,desc_data,info,work=aux,trans=trans_) + if (info /= psb_success_) exit - if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + call sm%sv%apply(cone,tx,cone,ty,desc_data,trans_,aux,wv(3:),info,init='Y') - if (info /= psb_success_) then - info=psb_err_internal_error_ - call psb_errpush(info,name,& - & a_err='subsolve with Jacobi sweeps > 1') - goto 9999 + if (info /= psb_success_) exit + + if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then + call psb_geaxpby(cone,x,czero,r,r,desc_data,info) + call psb_spmm(-cone,sm%pa,ty,cone,r,desc_data,info) + res = psb_genrm2(r,desc_data,info) + if( sm%printres ) then + call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol) + end if + if ( res < sm%tol*resdenum ) then + if( (sm%printres).and.(mod(sm%printiter,sm%checkiter)/=0) ) & + & call log_conv("BJAC",me,i,1,res,resdenum,sm%tol) + exit + end if end if + end do - end associate - class default - ! - ! - ! Apply multiple sweeps of a block-Jacobi solver - ! to compute an approximate solution of a linear system. - ! - ! - if (size(wv) < 2) then - info = psb_err_internal_error_ + if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + + if (info /= psb_success_) then + info=psb_err_internal_error_ call psb_errpush(info,name,& - & a_err='invalid wv size in smoother_apply') + & a_err='subsolve with Jacobi sweeps > 1') goto 9999 end if - associate(tx => wv(1), ty => wv(2)) - - ! - ! Unroll the first iteration and fold it inside SELECT CASE - ! this will save one AXPBY and one SPMM when INIT=Z, and will be - ! significant when sweeps=1 (a common case) - ! - select case (init_) - case('Z') - - call sm%sv%apply(cone,x,czero,ty,desc_data,trans_,aux,wv(3:),info,init='Z') - - case('Y') - call psb_geaxpby(cone,x,czero,tx,desc_data,info) - call psb_geaxpby(cone,y,czero,ty,desc_data,info) - call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,wv(3:),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,x,czero,tx,desc_data,info) - call psb_geaxpby(cone,initu,czero,ty,desc_data,info) - call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='wrong init to smoother_apply') - goto 9999 - end select - - 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. - ! - call psb_geaxpby(cone,x,czero,tx,desc_data,info) - call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') - - if (info /= psb_success_) exit - - if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then - call psb_geaxpby(cone,x,czero,r,r,desc_data,info) - call psb_spmm(-cone,sm%pa,ty,cone,r,desc_data,info) - res = psb_genrm2(r,desc_data,info) - if( sm%printres ) then - call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol) - end if - if (res < sm%tol*resdenum ) then - if( (sm%printres).and.( mod(sm%printiter,sm%checkiter) /=0 ) ) & - & call log_conv("BJAC",me,i,1,res,resdenum,sm%tol) - exit - end if - end if - - end do - - if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) - - 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 + end associate - end associate - end select + else if (sweeps == 0) then + ! + ! 0 sweeps of smoother is the identity operator + ! + call psb_geaxpby(alpha,x,beta,y,desc_data,info) else diff --git a/mlprec/impl/smoother/mld_c_jac_smoother_bld.f90 b/mlprec/impl/smoother/mld_c_jac_smoother_bld.f90 index c443886a..81c1feb8 100644 --- a/mlprec/impl/smoother/mld_c_jac_smoother_bld.f90 +++ b/mlprec/impl/smoother/mld_c_jac_smoother_bld.f90 @@ -71,13 +71,12 @@ subroutine mld_c_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold) nrow_a = a%get_nrows() nztota = a%get_nzeros() - if( sm%checkres ) sm%pa => a + sm%pa => a select type (smsv => sm%sv) class is (mld_c_diag_solver_type) call sm%nd%free() - sm%pa => a - sm%nd_nnz_tot = nztota + sm%nd_nnz_tot = nztota call psb_sum(ictxt,sm%nd_nnz_tot) call sm%sv%build(a,desc_a,info,amold=amold,vmold=vmold) diff --git a/mlprec/impl/smoother/mld_c_jac_smoother_clone.f90 b/mlprec/impl/smoother/mld_c_jac_smoother_clone.f90 index f17f50a6..9b52cc60 100644 --- a/mlprec/impl/smoother/mld_c_jac_smoother_clone.f90 +++ b/mlprec/impl/smoother/mld_c_jac_smoother_clone.f90 @@ -77,6 +77,7 @@ subroutine mld_c_jac_smoother_clone(sm,smout,info) allocate(smout%sv,mold=sm%sv,stat=info) if (info == psb_success_) call sm%sv%clone(smo%sv,info) end if + smo%pa => sm%pa class default info = psb_err_internal_error_ diff --git a/mlprec/impl/smoother/mld_d_as_smoother_bld.f90 b/mlprec/impl/smoother/mld_d_as_smoother_bld.f90 index 6cc55934..56497eae 100644 --- a/mlprec/impl/smoother/mld_d_as_smoother_bld.f90 +++ b/mlprec/impl/smoother/mld_d_as_smoother_bld.f90 @@ -66,7 +66,7 @@ subroutine mld_d_as_smoother_bld(a,desc_a,sm,info,amold,vmold,imold) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),' start' - + sm%pa => a novr = sm%novr if (novr < 0) then info=psb_err_invalid_ovr_num_ diff --git a/mlprec/impl/smoother/mld_d_as_smoother_clone.f90 b/mlprec/impl/smoother/mld_d_as_smoother_clone.f90 index e84e670c..7703a6fd 100644 --- a/mlprec/impl/smoother/mld_d_as_smoother_clone.f90 +++ b/mlprec/impl/smoother/mld_d_as_smoother_clone.f90 @@ -77,6 +77,7 @@ subroutine mld_d_as_smoother_clone(sm,smout,info) allocate(smout%sv,mold=sm%sv,stat=info) if (info == psb_success_) call sm%sv%clone(smo%sv,info) end if + smo%pa => sm%pa class default info = psb_err_internal_error_ diff --git a/mlprec/impl/smoother/mld_d_as_smoother_free.f90 b/mlprec/impl/smoother/mld_d_as_smoother_free.f90 index 2b446d9d..1a56ecb8 100644 --- a/mlprec/impl/smoother/mld_d_as_smoother_free.f90 +++ b/mlprec/impl/smoother/mld_d_as_smoother_free.f90 @@ -61,6 +61,7 @@ subroutine mld_d_as_smoother_free(sm,info) end if end if call sm%nd%free() + sm%pa => null() call psb_erractionrestore(err_act) return diff --git a/mlprec/impl/smoother/mld_d_jac_smoother_apply.f90 b/mlprec/impl/smoother/mld_d_jac_smoother_apply.f90 index a2512551..651dc892 100644 --- a/mlprec/impl/smoother/mld_d_jac_smoother_apply.f90 +++ b/mlprec/impl/smoother/mld_d_jac_smoother_apply.f90 @@ -65,7 +65,7 @@ subroutine mld_d_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,& ictxt = desc_data%get_context() call psb_info(ictxt,me,np) - + if (present(init)) then init_ = psb_toupper(init) else @@ -102,146 +102,63 @@ subroutine mld_d_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,& goto 9999 end if endif - - if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then - ! if .not.sv%is_iterative, there's no need to pass init - call sm%sv%apply(alpha,x,beta,y,desc_data,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 - - else if (sweeps >= 0) then - if (associated(sm%pa)) then - ! - ! This means we are dealing with a pure Jacobi smoother/solver. - ! - call psb_geasb(tx,desc_data,info) - call psb_geasb(ty,desc_data,info) - select case (init_) - case('Z') - - call sm%sv%apply(done,x,dzero,ty,desc_data,trans_,aux,info,init='Z') - - case('Y') - call psb_geaxpby(done,x,dzero,tx,desc_data,info) - call psb_geaxpby(done,y,dzero,ty,desc_data,info) - call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(done,tx,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,x,dzero,tx,desc_data,info) - call psb_geaxpby(done,initu,dzero,ty,desc_data,info) - call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,info,init='Y') - - case default + + if (sweeps > 0) then + call psb_geasb(tx,desc_data,info) + call psb_geasb(ty,desc_data,info) + select case (init_) + case('Z') + + call sm%sv%apply(done,x,dzero,ty,desc_data,trans_,aux,info,init='Z') + + case('Y') + call psb_geaxpby(done,x,dzero,tx,desc_data,info) + call psb_geaxpby(done,y,dzero,ty,desc_data,info) + call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(done,tx,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='wrong init to smoother_apply') + & a_err='missing initu to smoother_apply') goto 9999 - end select - - do i=1, sweeps-1 - ! - ! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)), - ! where is the diagonal and A the matrix. - ! - call psb_geaxpby(done,x,dzero,tx,desc_data,info) - call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call sm%sv%apply(done,tx,done,ty,desc_data,trans_,aux,info,init='Y') - - if (info /= psb_success_) exit - end do - - - if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) - - 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 - ! - ! - ! Apply multiple sweeps of a block-Jacobi solver - ! to compute an approximate solution of a linear system. - ! - ! - call psb_geasb(tx,desc_data,info) - call psb_geasb(ty,desc_data,info) + call psb_geaxpby(done,x,dzero,tx,desc_data,info) + call psb_geaxpby(done,initu,dzero,ty,desc_data,info) + call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(done,tx,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 + do i=1, sweeps-1 ! - ! Unroll the first iteration and fold it inside SELECT CASE - ! this will save one AXPBY and one SPMM when INIT=Z, and will be - ! significant when sweeps=1 (a common case) + ! Compute Y(j+1) = Y(j)+ M^(-1)*(X-A*Y(j)), ! - select case (init_) - case('Z') - - call sm%sv%apply(done,x,dzero,ty,desc_data,trans_,aux,info,init='Z') - - case('Y') - call psb_geaxpby(done,x,dzero,tx,desc_data,info) - call psb_geaxpby(done,y,dzero,ty,desc_data,info) - call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(done,tx,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,x,dzero,tx,desc_data,info) - call psb_geaxpby(done,initu,dzero,ty,desc_data,info) - call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(done,tx,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 + call psb_geaxpby(done,x,dzero,tx,desc_data,info) + call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_) - 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. - ! - call psb_geaxpby(done,x,dzero,tx,desc_data,info) - call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_) + if (info /= psb_success_) exit - if (info /= psb_success_) exit + call sm%sv%apply(done,tx,done,ty,desc_data,trans_,aux,info,init='Y') - call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,info,init='Y') + if (info /= psb_success_) exit + end do - if (info /= psb_success_) exit - end do - if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) - 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 + 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 - + + deallocate(tx,ty,stat=info) if (info /= psb_success_) then info=psb_err_internal_error_ @@ -250,6 +167,12 @@ subroutine mld_d_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,& goto 9999 end if + else if (sweeps == 0) then + ! + ! 0 sweeps of smoother is the identity operator + ! + call psb_geaxpby(alpha,x,beta,y,desc_data,info) + else info = psb_err_iarg_neg_ diff --git a/mlprec/impl/smoother/mld_d_jac_smoother_apply_vect.f90 b/mlprec/impl/smoother/mld_d_jac_smoother_apply_vect.f90 index 2e95454d..c243d672 100644 --- a/mlprec/impl/smoother/mld_d_jac_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_d_jac_smoother_apply_vect.f90 @@ -49,7 +49,7 @@ subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& type(psb_d_vect_type),intent(inout) :: y real(psb_dpk_),intent(in) :: alpha,beta character(len=1),intent(in) :: trans - integer(psb_ipk_), intent(in) :: sweeps + integer(psb_ipk_), intent(in) :: sweeps real(psb_dpk_),target, intent(inout) :: work(:) type(psb_d_vect_type),intent(inout) :: wv(:) integer(psb_ipk_), intent(out) :: info @@ -108,190 +108,90 @@ subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& end if endif - if(sm%checkres) then - call psb_geall(r,desc_data,info) - call psb_geasb(r,desc_data,info) - resdenum = psb_genrm2(x,desc_data,info) - end if + if (sweeps > 0) then - if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then - ! if .not.sv%is_iterative, there's no need to pass init - call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info) + if(sm%checkres) then + call psb_geall(r,desc_data,info) + call psb_geasb(r,desc_data,info) + resdenum = psb_genrm2(x,desc_data,info) + end if - 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 - - else if (sweeps >= 0) then - select type (smsv => sm%sv) - class is (mld_d_diag_solver_type) - ! - ! This means we are dealing with a pure Jacobi smoother/solver. - ! - associate(tx => wv(1), ty => wv(2)) - select case (init_) - case('Z') - - call sm%sv%apply(done,x,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Z') - - case('Y') - call psb_geaxpby(done,x,dzero,tx,desc_data,info) - call psb_geaxpby(done,y,dzero,ty,desc_data,info) - call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,wv(3:),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,x,dzero,tx,desc_data,info) - call psb_geaxpby(done,initu,dzero,ty,desc_data,info) - call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') + associate(tx => wv(1), ty => wv(2)) + select case (init_) + case('Z') + + call sm%sv%apply(done,x,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Z') + + case('Y') + call psb_geaxpby(done,x,dzero,tx,desc_data,info) + call psb_geaxpby(done,y,dzero,ty,desc_data,info) + call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') - case default + case('U') + if (.not.present(initu)) then call psb_errpush(psb_err_internal_error_,name,& - & a_err='wrong init to smoother_apply') + & a_err='missing initu to smoother_apply') goto 9999 - end select - - do i=1, sweeps-1 - ! - ! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)), - ! where is the diagonal and A the matrix. - ! - call psb_geaxpby(done,x,dzero,tx,desc_data,info) - call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call sm%sv%apply(done,tx,done,ty,desc_data,trans_,aux,wv(3:),info,init='Y') - - if (info /= psb_success_) exit - - if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then - call psb_geaxpby(done,x,dzero,r,r,desc_data,info) - call psb_spmm(-done,sm%pa,ty,done,r,desc_data,info) - res = psb_genrm2(r,desc_data,info) - if( sm%printres ) then - call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol) - end if - if ( res < sm%tol*resdenum ) then - if( (sm%printres).and.(mod(sm%printiter,sm%checkiter)/=0) ) & - & call log_conv("BJAC",me,i,1,res,resdenum,sm%tol) - exit - end if - end if + end if + call psb_geaxpby(done,x,dzero,tx,desc_data,info) + call psb_geaxpby(done,initu,dzero,ty,desc_data,info) + call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') + + case default + call psb_errpush(psb_err_internal_error_,name,& + & a_err='wrong init to smoother_apply') + goto 9999 + end select - end do + do i=1, sweeps-1 + ! + ! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)), + ! where is the diagonal and A the matrix. + ! + call psb_geaxpby(done,x,dzero,tx,desc_data,info) + call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_) + if (info /= psb_success_) exit - if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + call sm%sv%apply(done,tx,done,ty,desc_data,trans_,aux,wv(3:),info,init='Y') - if (info /= psb_success_) then - info=psb_err_internal_error_ - call psb_errpush(info,name,& - & a_err='subsolve with Jacobi sweeps > 1') - goto 9999 + if (info /= psb_success_) exit + + if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then + call psb_geaxpby(done,x,dzero,r,r,desc_data,info) + call psb_spmm(-done,sm%pa,ty,done,r,desc_data,info) + res = psb_genrm2(r,desc_data,info) + if( sm%printres ) then + call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol) + end if + if ( res < sm%tol*resdenum ) then + if( (sm%printres).and.(mod(sm%printiter,sm%checkiter)/=0) ) & + & call log_conv("BJAC",me,i,1,res,resdenum,sm%tol) + exit + end if end if + end do - end associate - class default - ! - ! - ! Apply multiple sweeps of a block-Jacobi solver - ! to compute an approximate solution of a linear system. - ! - ! - if (size(wv) < 2) then - info = psb_err_internal_error_ + if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + + if (info /= psb_success_) then + info=psb_err_internal_error_ call psb_errpush(info,name,& - & a_err='invalid wv size in smoother_apply') + & a_err='subsolve with Jacobi sweeps > 1') goto 9999 end if - associate(tx => wv(1), ty => wv(2)) - - ! - ! Unroll the first iteration and fold it inside SELECT CASE - ! this will save one AXPBY and one SPMM when INIT=Z, and will be - ! significant when sweeps=1 (a common case) - ! - select case (init_) - case('Z') - - call sm%sv%apply(done,x,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Z') - - case('Y') - call psb_geaxpby(done,x,dzero,tx,desc_data,info) - call psb_geaxpby(done,y,dzero,ty,desc_data,info) - call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,wv(3:),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,x,dzero,tx,desc_data,info) - call psb_geaxpby(done,initu,dzero,ty,desc_data,info) - call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='wrong init to smoother_apply') - goto 9999 - end select - - 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. - ! - call psb_geaxpby(done,x,dzero,tx,desc_data,info) - call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') - - if (info /= psb_success_) exit - - if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then - call psb_geaxpby(done,x,dzero,r,r,desc_data,info) - call psb_spmm(-done,sm%pa,ty,done,r,desc_data,info) - res = psb_genrm2(r,desc_data,info) - if( sm%printres ) then - call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol) - end if - if (res < sm%tol*resdenum ) then - if( (sm%printres).and.( mod(sm%printiter,sm%checkiter) /=0 ) ) & - & call log_conv("BJAC",me,i,1,res,resdenum,sm%tol) - exit - end if - end if - - end do - - if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) - - 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 + end associate - end associate - end select + else if (sweeps == 0) then + ! + ! 0 sweeps of smoother is the identity operator + ! + call psb_geaxpby(alpha,x,beta,y,desc_data,info) else diff --git a/mlprec/impl/smoother/mld_d_jac_smoother_bld.f90 b/mlprec/impl/smoother/mld_d_jac_smoother_bld.f90 index b008957a..b84396a6 100644 --- a/mlprec/impl/smoother/mld_d_jac_smoother_bld.f90 +++ b/mlprec/impl/smoother/mld_d_jac_smoother_bld.f90 @@ -71,13 +71,12 @@ subroutine mld_d_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold) nrow_a = a%get_nrows() nztota = a%get_nzeros() - if( sm%checkres ) sm%pa => a + sm%pa => a select type (smsv => sm%sv) class is (mld_d_diag_solver_type) call sm%nd%free() - sm%pa => a - sm%nd_nnz_tot = nztota + sm%nd_nnz_tot = nztota call psb_sum(ictxt,sm%nd_nnz_tot) call sm%sv%build(a,desc_a,info,amold=amold,vmold=vmold) diff --git a/mlprec/impl/smoother/mld_d_jac_smoother_clone.f90 b/mlprec/impl/smoother/mld_d_jac_smoother_clone.f90 index caa88534..e4b8b77e 100644 --- a/mlprec/impl/smoother/mld_d_jac_smoother_clone.f90 +++ b/mlprec/impl/smoother/mld_d_jac_smoother_clone.f90 @@ -77,6 +77,7 @@ subroutine mld_d_jac_smoother_clone(sm,smout,info) allocate(smout%sv,mold=sm%sv,stat=info) if (info == psb_success_) call sm%sv%clone(smo%sv,info) end if + smo%pa => sm%pa class default info = psb_err_internal_error_ diff --git a/mlprec/impl/smoother/mld_s_as_smoother_bld.f90 b/mlprec/impl/smoother/mld_s_as_smoother_bld.f90 index f5f70455..ab13d02a 100644 --- a/mlprec/impl/smoother/mld_s_as_smoother_bld.f90 +++ b/mlprec/impl/smoother/mld_s_as_smoother_bld.f90 @@ -66,7 +66,7 @@ subroutine mld_s_as_smoother_bld(a,desc_a,sm,info,amold,vmold,imold) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),' start' - + sm%pa => a novr = sm%novr if (novr < 0) then info=psb_err_invalid_ovr_num_ diff --git a/mlprec/impl/smoother/mld_s_as_smoother_clone.f90 b/mlprec/impl/smoother/mld_s_as_smoother_clone.f90 index 6076c287..2dc9d2cf 100644 --- a/mlprec/impl/smoother/mld_s_as_smoother_clone.f90 +++ b/mlprec/impl/smoother/mld_s_as_smoother_clone.f90 @@ -77,6 +77,7 @@ subroutine mld_s_as_smoother_clone(sm,smout,info) allocate(smout%sv,mold=sm%sv,stat=info) if (info == psb_success_) call sm%sv%clone(smo%sv,info) end if + smo%pa => sm%pa class default info = psb_err_internal_error_ diff --git a/mlprec/impl/smoother/mld_s_as_smoother_free.f90 b/mlprec/impl/smoother/mld_s_as_smoother_free.f90 index 5a632a42..cc34b817 100644 --- a/mlprec/impl/smoother/mld_s_as_smoother_free.f90 +++ b/mlprec/impl/smoother/mld_s_as_smoother_free.f90 @@ -61,6 +61,7 @@ subroutine mld_s_as_smoother_free(sm,info) end if end if call sm%nd%free() + sm%pa => null() call psb_erractionrestore(err_act) return diff --git a/mlprec/impl/smoother/mld_s_jac_smoother_apply.f90 b/mlprec/impl/smoother/mld_s_jac_smoother_apply.f90 index a7b28c29..06cf931e 100644 --- a/mlprec/impl/smoother/mld_s_jac_smoother_apply.f90 +++ b/mlprec/impl/smoother/mld_s_jac_smoother_apply.f90 @@ -65,7 +65,7 @@ subroutine mld_s_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,& ictxt = desc_data%get_context() call psb_info(ictxt,me,np) - + if (present(init)) then init_ = psb_toupper(init) else @@ -102,146 +102,63 @@ subroutine mld_s_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,& goto 9999 end if endif - - if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then - ! if .not.sv%is_iterative, there's no need to pass init - call sm%sv%apply(alpha,x,beta,y,desc_data,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 - - else if (sweeps >= 0) then - if (associated(sm%pa)) then - ! - ! This means we are dealing with a pure Jacobi smoother/solver. - ! - call psb_geasb(tx,desc_data,info) - call psb_geasb(ty,desc_data,info) - select case (init_) - case('Z') - - call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,info,init='Z') - - case('Y') - call psb_geaxpby(sone,x,szero,tx,desc_data,info) - call psb_geaxpby(sone,y,szero,ty,desc_data,info) - call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(sone,tx,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,x,szero,tx,desc_data,info) - call psb_geaxpby(sone,initu,szero,ty,desc_data,info) - call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,info,init='Y') - - case default + + if (sweeps > 0) then + call psb_geasb(tx,desc_data,info) + call psb_geasb(ty,desc_data,info) + select case (init_) + case('Z') + + call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,info,init='Z') + + case('Y') + call psb_geaxpby(sone,x,szero,tx,desc_data,info) + call psb_geaxpby(sone,y,szero,ty,desc_data,info) + call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(sone,tx,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='wrong init to smoother_apply') + & a_err='missing initu to smoother_apply') goto 9999 - end select - - do i=1, sweeps-1 - ! - ! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)), - ! where is the diagonal and A the matrix. - ! - call psb_geaxpby(sone,x,szero,tx,desc_data,info) - call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call sm%sv%apply(sone,tx,sone,ty,desc_data,trans_,aux,info,init='Y') - - if (info /= psb_success_) exit - end do - - - if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) - - 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 - ! - ! - ! Apply multiple sweeps of a block-Jacobi solver - ! to compute an approximate solution of a linear system. - ! - ! - call psb_geasb(tx,desc_data,info) - call psb_geasb(ty,desc_data,info) + call psb_geaxpby(sone,x,szero,tx,desc_data,info) + call psb_geaxpby(sone,initu,szero,ty,desc_data,info) + call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(sone,tx,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 + do i=1, sweeps-1 ! - ! Unroll the first iteration and fold it inside SELECT CASE - ! this will save one AXPBY and one SPMM when INIT=Z, and will be - ! significant when sweeps=1 (a common case) + ! Compute Y(j+1) = Y(j)+ M^(-1)*(X-A*Y(j)), ! - select case (init_) - case('Z') - - call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,info,init='Z') - - case('Y') - call psb_geaxpby(sone,x,szero,tx,desc_data,info) - call psb_geaxpby(sone,y,szero,ty,desc_data,info) - call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(sone,tx,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,x,szero,tx,desc_data,info) - call psb_geaxpby(sone,initu,szero,ty,desc_data,info) - call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(sone,tx,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 + call psb_geaxpby(sone,x,szero,tx,desc_data,info) + call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_) - 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. - ! - call psb_geaxpby(sone,x,szero,tx,desc_data,info) - call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_) + if (info /= psb_success_) exit - if (info /= psb_success_) exit + call sm%sv%apply(sone,tx,sone,ty,desc_data,trans_,aux,info,init='Y') - call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,info,init='Y') + if (info /= psb_success_) exit + end do - if (info /= psb_success_) exit - end do - if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) - 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 + 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 - + + deallocate(tx,ty,stat=info) if (info /= psb_success_) then info=psb_err_internal_error_ @@ -250,6 +167,12 @@ subroutine mld_s_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,& goto 9999 end if + else if (sweeps == 0) then + ! + ! 0 sweeps of smoother is the identity operator + ! + call psb_geaxpby(alpha,x,beta,y,desc_data,info) + else info = psb_err_iarg_neg_ diff --git a/mlprec/impl/smoother/mld_s_jac_smoother_apply_vect.f90 b/mlprec/impl/smoother/mld_s_jac_smoother_apply_vect.f90 index 5ed28511..6560cfa9 100644 --- a/mlprec/impl/smoother/mld_s_jac_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_s_jac_smoother_apply_vect.f90 @@ -49,7 +49,7 @@ subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& type(psb_s_vect_type),intent(inout) :: y real(psb_spk_),intent(in) :: alpha,beta character(len=1),intent(in) :: trans - integer(psb_ipk_), intent(in) :: sweeps + integer(psb_ipk_), intent(in) :: sweeps real(psb_spk_),target, intent(inout) :: work(:) type(psb_s_vect_type),intent(inout) :: wv(:) integer(psb_ipk_), intent(out) :: info @@ -108,190 +108,90 @@ subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& end if endif - if(sm%checkres) then - call psb_geall(r,desc_data,info) - call psb_geasb(r,desc_data,info) - resdenum = psb_genrm2(x,desc_data,info) - end if + if (sweeps > 0) then - if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then - ! if .not.sv%is_iterative, there's no need to pass init - call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info) + if(sm%checkres) then + call psb_geall(r,desc_data,info) + call psb_geasb(r,desc_data,info) + resdenum = psb_genrm2(x,desc_data,info) + end if - 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 - - else if (sweeps >= 0) then - select type (smsv => sm%sv) - class is (mld_s_diag_solver_type) - ! - ! This means we are dealing with a pure Jacobi smoother/solver. - ! - associate(tx => wv(1), ty => wv(2)) - select case (init_) - case('Z') - - call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Z') - - case('Y') - call psb_geaxpby(sone,x,szero,tx,desc_data,info) - call psb_geaxpby(sone,y,szero,ty,desc_data,info) - call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),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,x,szero,tx,desc_data,info) - call psb_geaxpby(sone,initu,szero,ty,desc_data,info) - call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') + associate(tx => wv(1), ty => wv(2)) + select case (init_) + case('Z') + + call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Z') + + case('Y') + call psb_geaxpby(sone,x,szero,tx,desc_data,info) + call psb_geaxpby(sone,y,szero,ty,desc_data,info) + call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') - case default + case('U') + if (.not.present(initu)) then call psb_errpush(psb_err_internal_error_,name,& - & a_err='wrong init to smoother_apply') + & a_err='missing initu to smoother_apply') goto 9999 - end select - - do i=1, sweeps-1 - ! - ! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)), - ! where is the diagonal and A the matrix. - ! - call psb_geaxpby(sone,x,szero,tx,desc_data,info) - call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call sm%sv%apply(sone,tx,sone,ty,desc_data,trans_,aux,wv(3:),info,init='Y') - - if (info /= psb_success_) exit - - if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then - call psb_geaxpby(sone,x,szero,r,r,desc_data,info) - call psb_spmm(-sone,sm%pa,ty,sone,r,desc_data,info) - res = psb_genrm2(r,desc_data,info) - if( sm%printres ) then - call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol) - end if - if ( res < sm%tol*resdenum ) then - if( (sm%printres).and.(mod(sm%printiter,sm%checkiter)/=0) ) & - & call log_conv("BJAC",me,i,1,res,resdenum,sm%tol) - exit - end if - end if + end if + call psb_geaxpby(sone,x,szero,tx,desc_data,info) + call psb_geaxpby(sone,initu,szero,ty,desc_data,info) + call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') + + case default + call psb_errpush(psb_err_internal_error_,name,& + & a_err='wrong init to smoother_apply') + goto 9999 + end select - end do + do i=1, sweeps-1 + ! + ! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)), + ! where is the diagonal and A the matrix. + ! + call psb_geaxpby(sone,x,szero,tx,desc_data,info) + call psb_spmm(-sone,sm%pa,ty,sone,tx,desc_data,info,work=aux,trans=trans_) + if (info /= psb_success_) exit - if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + call sm%sv%apply(sone,tx,sone,ty,desc_data,trans_,aux,wv(3:),info,init='Y') - if (info /= psb_success_) then - info=psb_err_internal_error_ - call psb_errpush(info,name,& - & a_err='subsolve with Jacobi sweeps > 1') - goto 9999 + if (info /= psb_success_) exit + + if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then + call psb_geaxpby(sone,x,szero,r,r,desc_data,info) + call psb_spmm(-sone,sm%pa,ty,sone,r,desc_data,info) + res = psb_genrm2(r,desc_data,info) + if( sm%printres ) then + call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol) + end if + if ( res < sm%tol*resdenum ) then + if( (sm%printres).and.(mod(sm%printiter,sm%checkiter)/=0) ) & + & call log_conv("BJAC",me,i,1,res,resdenum,sm%tol) + exit + end if end if + end do - end associate - class default - ! - ! - ! Apply multiple sweeps of a block-Jacobi solver - ! to compute an approximate solution of a linear system. - ! - ! - if (size(wv) < 2) then - info = psb_err_internal_error_ + if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + + if (info /= psb_success_) then + info=psb_err_internal_error_ call psb_errpush(info,name,& - & a_err='invalid wv size in smoother_apply') + & a_err='subsolve with Jacobi sweeps > 1') goto 9999 end if - associate(tx => wv(1), ty => wv(2)) - - ! - ! Unroll the first iteration and fold it inside SELECT CASE - ! this will save one AXPBY and one SPMM when INIT=Z, and will be - ! significant when sweeps=1 (a common case) - ! - select case (init_) - case('Z') - - call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Z') - - case('Y') - call psb_geaxpby(sone,x,szero,tx,desc_data,info) - call psb_geaxpby(sone,y,szero,ty,desc_data,info) - call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),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,x,szero,tx,desc_data,info) - call psb_geaxpby(sone,initu,szero,ty,desc_data,info) - call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='wrong init to smoother_apply') - goto 9999 - end select - - 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. - ! - call psb_geaxpby(sone,x,szero,tx,desc_data,info) - call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') - - if (info /= psb_success_) exit - - if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then - call psb_geaxpby(sone,x,szero,r,r,desc_data,info) - call psb_spmm(-sone,sm%pa,ty,sone,r,desc_data,info) - res = psb_genrm2(r,desc_data,info) - if( sm%printres ) then - call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol) - end if - if (res < sm%tol*resdenum ) then - if( (sm%printres).and.( mod(sm%printiter,sm%checkiter) /=0 ) ) & - & call log_conv("BJAC",me,i,1,res,resdenum,sm%tol) - exit - end if - end if - - end do - - if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) - - 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 + end associate - end associate - end select + else if (sweeps == 0) then + ! + ! 0 sweeps of smoother is the identity operator + ! + call psb_geaxpby(alpha,x,beta,y,desc_data,info) else diff --git a/mlprec/impl/smoother/mld_s_jac_smoother_bld.f90 b/mlprec/impl/smoother/mld_s_jac_smoother_bld.f90 index deb7acae..deedaaa6 100644 --- a/mlprec/impl/smoother/mld_s_jac_smoother_bld.f90 +++ b/mlprec/impl/smoother/mld_s_jac_smoother_bld.f90 @@ -71,13 +71,12 @@ subroutine mld_s_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold) nrow_a = a%get_nrows() nztota = a%get_nzeros() - if( sm%checkres ) sm%pa => a + sm%pa => a select type (smsv => sm%sv) class is (mld_s_diag_solver_type) call sm%nd%free() - sm%pa => a - sm%nd_nnz_tot = nztota + sm%nd_nnz_tot = nztota call psb_sum(ictxt,sm%nd_nnz_tot) call sm%sv%build(a,desc_a,info,amold=amold,vmold=vmold) diff --git a/mlprec/impl/smoother/mld_s_jac_smoother_clone.f90 b/mlprec/impl/smoother/mld_s_jac_smoother_clone.f90 index a0b6c349..abfb7660 100644 --- a/mlprec/impl/smoother/mld_s_jac_smoother_clone.f90 +++ b/mlprec/impl/smoother/mld_s_jac_smoother_clone.f90 @@ -77,6 +77,7 @@ subroutine mld_s_jac_smoother_clone(sm,smout,info) allocate(smout%sv,mold=sm%sv,stat=info) if (info == psb_success_) call sm%sv%clone(smo%sv,info) end if + smo%pa => sm%pa class default info = psb_err_internal_error_ diff --git a/mlprec/impl/smoother/mld_z_as_smoother_bld.f90 b/mlprec/impl/smoother/mld_z_as_smoother_bld.f90 index b002c1af..703d3db5 100644 --- a/mlprec/impl/smoother/mld_z_as_smoother_bld.f90 +++ b/mlprec/impl/smoother/mld_z_as_smoother_bld.f90 @@ -66,7 +66,7 @@ subroutine mld_z_as_smoother_bld(a,desc_a,sm,info,amold,vmold,imold) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),' start' - + sm%pa => a novr = sm%novr if (novr < 0) then info=psb_err_invalid_ovr_num_ diff --git a/mlprec/impl/smoother/mld_z_as_smoother_clone.f90 b/mlprec/impl/smoother/mld_z_as_smoother_clone.f90 index c192cdc6..542fade9 100644 --- a/mlprec/impl/smoother/mld_z_as_smoother_clone.f90 +++ b/mlprec/impl/smoother/mld_z_as_smoother_clone.f90 @@ -77,6 +77,7 @@ subroutine mld_z_as_smoother_clone(sm,smout,info) allocate(smout%sv,mold=sm%sv,stat=info) if (info == psb_success_) call sm%sv%clone(smo%sv,info) end if + smo%pa => sm%pa class default info = psb_err_internal_error_ diff --git a/mlprec/impl/smoother/mld_z_as_smoother_free.f90 b/mlprec/impl/smoother/mld_z_as_smoother_free.f90 index 6745ff10..a279d80c 100644 --- a/mlprec/impl/smoother/mld_z_as_smoother_free.f90 +++ b/mlprec/impl/smoother/mld_z_as_smoother_free.f90 @@ -61,6 +61,7 @@ subroutine mld_z_as_smoother_free(sm,info) end if end if call sm%nd%free() + sm%pa => null() call psb_erractionrestore(err_act) return diff --git a/mlprec/impl/smoother/mld_z_jac_smoother_apply.f90 b/mlprec/impl/smoother/mld_z_jac_smoother_apply.f90 index 3b507f5f..6e13f032 100644 --- a/mlprec/impl/smoother/mld_z_jac_smoother_apply.f90 +++ b/mlprec/impl/smoother/mld_z_jac_smoother_apply.f90 @@ -65,7 +65,7 @@ subroutine mld_z_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,& ictxt = desc_data%get_context() call psb_info(ictxt,me,np) - + if (present(init)) then init_ = psb_toupper(init) else @@ -102,146 +102,63 @@ subroutine mld_z_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,& goto 9999 end if endif - - if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then - ! if .not.sv%is_iterative, there's no need to pass init - call sm%sv%apply(alpha,x,beta,y,desc_data,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 - - else if (sweeps >= 0) then - if (associated(sm%pa)) then - ! - ! This means we are dealing with a pure Jacobi smoother/solver. - ! - call psb_geasb(tx,desc_data,info) - call psb_geasb(ty,desc_data,info) - select case (init_) - case('Z') - - call sm%sv%apply(zone,x,zzero,ty,desc_data,trans_,aux,info,init='Z') - - case('Y') - call psb_geaxpby(zone,x,zzero,tx,desc_data,info) - call psb_geaxpby(zone,y,zzero,ty,desc_data,info) - call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(zone,tx,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,x,zzero,tx,desc_data,info) - call psb_geaxpby(zone,initu,zzero,ty,desc_data,info) - call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,info,init='Y') - - case default + + if (sweeps > 0) then + call psb_geasb(tx,desc_data,info) + call psb_geasb(ty,desc_data,info) + select case (init_) + case('Z') + + call sm%sv%apply(zone,x,zzero,ty,desc_data,trans_,aux,info,init='Z') + + case('Y') + call psb_geaxpby(zone,x,zzero,tx,desc_data,info) + call psb_geaxpby(zone,y,zzero,ty,desc_data,info) + call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(zone,tx,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='wrong init to smoother_apply') + & a_err='missing initu to smoother_apply') goto 9999 - end select - - do i=1, sweeps-1 - ! - ! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)), - ! where is the diagonal and A the matrix. - ! - call psb_geaxpby(zone,x,zzero,tx,desc_data,info) - call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call sm%sv%apply(zone,tx,zone,ty,desc_data,trans_,aux,info,init='Y') - - if (info /= psb_success_) exit - end do - - - if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) - - 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 - ! - ! - ! Apply multiple sweeps of a block-Jacobi solver - ! to compute an approximate solution of a linear system. - ! - ! - call psb_geasb(tx,desc_data,info) - call psb_geasb(ty,desc_data,info) + call psb_geaxpby(zone,x,zzero,tx,desc_data,info) + call psb_geaxpby(zone,initu,zzero,ty,desc_data,info) + call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(zone,tx,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 + do i=1, sweeps-1 ! - ! Unroll the first iteration and fold it inside SELECT CASE - ! this will save one AXPBY and one SPMM when INIT=Z, and will be - ! significant when sweeps=1 (a common case) + ! Compute Y(j+1) = Y(j)+ M^(-1)*(X-A*Y(j)), ! - select case (init_) - case('Z') - - call sm%sv%apply(zone,x,zzero,ty,desc_data,trans_,aux,info,init='Z') - - case('Y') - call psb_geaxpby(zone,x,zzero,tx,desc_data,info) - call psb_geaxpby(zone,y,zzero,ty,desc_data,info) - call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(zone,tx,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,x,zzero,tx,desc_data,info) - call psb_geaxpby(zone,initu,zzero,ty,desc_data,info) - call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(zone,tx,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 + call psb_geaxpby(zone,x,zzero,tx,desc_data,info) + call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_) - 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. - ! - call psb_geaxpby(zone,x,zzero,tx,desc_data,info) - call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_) + if (info /= psb_success_) exit - if (info /= psb_success_) exit + call sm%sv%apply(zone,tx,zone,ty,desc_data,trans_,aux,info,init='Y') - call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,info,init='Y') + if (info /= psb_success_) exit + end do - if (info /= psb_success_) exit - end do - if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) - 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 + 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 - + + deallocate(tx,ty,stat=info) if (info /= psb_success_) then info=psb_err_internal_error_ @@ -250,6 +167,12 @@ subroutine mld_z_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,& goto 9999 end if + else if (sweeps == 0) then + ! + ! 0 sweeps of smoother is the identity operator + ! + call psb_geaxpby(alpha,x,beta,y,desc_data,info) + else info = psb_err_iarg_neg_ diff --git a/mlprec/impl/smoother/mld_z_jac_smoother_apply_vect.f90 b/mlprec/impl/smoother/mld_z_jac_smoother_apply_vect.f90 index 2da1aa53..e2912103 100644 --- a/mlprec/impl/smoother/mld_z_jac_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_z_jac_smoother_apply_vect.f90 @@ -49,7 +49,7 @@ subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& type(psb_z_vect_type),intent(inout) :: y complex(psb_dpk_),intent(in) :: alpha,beta character(len=1),intent(in) :: trans - integer(psb_ipk_), intent(in) :: sweeps + integer(psb_ipk_), intent(in) :: sweeps complex(psb_dpk_),target, intent(inout) :: work(:) type(psb_z_vect_type),intent(inout) :: wv(:) integer(psb_ipk_), intent(out) :: info @@ -108,190 +108,90 @@ subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& end if endif - if(sm%checkres) then - call psb_geall(r,desc_data,info) - call psb_geasb(r,desc_data,info) - resdenum = psb_genrm2(x,desc_data,info) - end if + if (sweeps > 0) then - if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then - ! if .not.sv%is_iterative, there's no need to pass init - call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info) + if(sm%checkres) then + call psb_geall(r,desc_data,info) + call psb_geasb(r,desc_data,info) + resdenum = psb_genrm2(x,desc_data,info) + end if - 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 - - else if (sweeps >= 0) then - select type (smsv => sm%sv) - class is (mld_z_diag_solver_type) - ! - ! This means we are dealing with a pure Jacobi smoother/solver. - ! - associate(tx => wv(1), ty => wv(2)) - select case (init_) - case('Z') - - call sm%sv%apply(zone,x,zzero,ty,desc_data,trans_,aux,wv(3:),info,init='Z') - - case('Y') - call psb_geaxpby(zone,x,zzero,tx,desc_data,info) - call psb_geaxpby(zone,y,zzero,ty,desc_data,info) - call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,wv(3:),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,x,zzero,tx,desc_data,info) - call psb_geaxpby(zone,initu,zzero,ty,desc_data,info) - call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') + associate(tx => wv(1), ty => wv(2)) + select case (init_) + case('Z') + + call sm%sv%apply(zone,x,zzero,ty,desc_data,trans_,aux,wv(3:),info,init='Z') + + case('Y') + call psb_geaxpby(zone,x,zzero,tx,desc_data,info) + call psb_geaxpby(zone,y,zzero,ty,desc_data,info) + call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') - case default + case('U') + if (.not.present(initu)) then call psb_errpush(psb_err_internal_error_,name,& - & a_err='wrong init to smoother_apply') + & a_err='missing initu to smoother_apply') goto 9999 - end select - - do i=1, sweeps-1 - ! - ! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)), - ! where is the diagonal and A the matrix. - ! - call psb_geaxpby(zone,x,zzero,tx,desc_data,info) - call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call sm%sv%apply(zone,tx,zone,ty,desc_data,trans_,aux,wv(3:),info,init='Y') - - if (info /= psb_success_) exit - - if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then - call psb_geaxpby(zone,x,zzero,r,r,desc_data,info) - call psb_spmm(-zone,sm%pa,ty,zone,r,desc_data,info) - res = psb_genrm2(r,desc_data,info) - if( sm%printres ) then - call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol) - end if - if ( res < sm%tol*resdenum ) then - if( (sm%printres).and.(mod(sm%printiter,sm%checkiter)/=0) ) & - & call log_conv("BJAC",me,i,1,res,resdenum,sm%tol) - exit - end if - end if + end if + call psb_geaxpby(zone,x,zzero,tx,desc_data,info) + call psb_geaxpby(zone,initu,zzero,ty,desc_data,info) + call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_) + call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') + + case default + call psb_errpush(psb_err_internal_error_,name,& + & a_err='wrong init to smoother_apply') + goto 9999 + end select - end do + do i=1, sweeps-1 + ! + ! Compute Y(j+1) = Y(j)+ D^(-1)*(X-A*Y(j)), + ! where is the diagonal and A the matrix. + ! + call psb_geaxpby(zone,x,zzero,tx,desc_data,info) + call psb_spmm(-zone,sm%pa,ty,zone,tx,desc_data,info,work=aux,trans=trans_) + if (info /= psb_success_) exit - if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + call sm%sv%apply(zone,tx,zone,ty,desc_data,trans_,aux,wv(3:),info,init='Y') - if (info /= psb_success_) then - info=psb_err_internal_error_ - call psb_errpush(info,name,& - & a_err='subsolve with Jacobi sweeps > 1') - goto 9999 + if (info /= psb_success_) exit + + if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then + call psb_geaxpby(zone,x,zzero,r,r,desc_data,info) + call psb_spmm(-zone,sm%pa,ty,zone,r,desc_data,info) + res = psb_genrm2(r,desc_data,info) + if( sm%printres ) then + call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol) + end if + if ( res < sm%tol*resdenum ) then + if( (sm%printres).and.(mod(sm%printiter,sm%checkiter)/=0) ) & + & call log_conv("BJAC",me,i,1,res,resdenum,sm%tol) + exit + end if end if + end do - end associate - class default - ! - ! - ! Apply multiple sweeps of a block-Jacobi solver - ! to compute an approximate solution of a linear system. - ! - ! - if (size(wv) < 2) then - info = psb_err_internal_error_ + if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + + if (info /= psb_success_) then + info=psb_err_internal_error_ call psb_errpush(info,name,& - & a_err='invalid wv size in smoother_apply') + & a_err='subsolve with Jacobi sweeps > 1') goto 9999 end if - associate(tx => wv(1), ty => wv(2)) - - ! - ! Unroll the first iteration and fold it inside SELECT CASE - ! this will save one AXPBY and one SPMM when INIT=Z, and will be - ! significant when sweeps=1 (a common case) - ! - select case (init_) - case('Z') - - call sm%sv%apply(zone,x,zzero,ty,desc_data,trans_,aux,wv(3:),info,init='Z') - - case('Y') - call psb_geaxpby(zone,x,zzero,tx,desc_data,info) - call psb_geaxpby(zone,y,zzero,ty,desc_data,info) - call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,wv(3:),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,x,zzero,tx,desc_data,info) - call psb_geaxpby(zone,initu,zzero,ty,desc_data,info) - call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_) - call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') - case default - call psb_errpush(psb_err_internal_error_,name,& - & a_err='wrong init to smoother_apply') - goto 9999 - end select - - 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. - ! - call psb_geaxpby(zone,x,zzero,tx,desc_data,info) - call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_) - - if (info /= psb_success_) exit - - call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,wv(3:),info,init='Y') - - if (info /= psb_success_) exit - - if ( sm%checkres.and.(mod(i,sm%checkiter) == 0) ) then - call psb_geaxpby(zone,x,zzero,r,r,desc_data,info) - call psb_spmm(-zone,sm%pa,ty,zone,r,desc_data,info) - res = psb_genrm2(r,desc_data,info) - if( sm%printres ) then - call log_conv("BJAC",me,i,sm%printiter,res,resdenum,sm%tol) - end if - if (res < sm%tol*resdenum ) then - if( (sm%printres).and.( mod(sm%printiter,sm%checkiter) /=0 ) ) & - & call log_conv("BJAC",me,i,1,res,resdenum,sm%tol) - exit - end if - end if - - end do - - if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) - - 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 + end associate - end associate - end select + else if (sweeps == 0) then + ! + ! 0 sweeps of smoother is the identity operator + ! + call psb_geaxpby(alpha,x,beta,y,desc_data,info) else diff --git a/mlprec/impl/smoother/mld_z_jac_smoother_bld.f90 b/mlprec/impl/smoother/mld_z_jac_smoother_bld.f90 index 3852f032..b52b7b62 100644 --- a/mlprec/impl/smoother/mld_z_jac_smoother_bld.f90 +++ b/mlprec/impl/smoother/mld_z_jac_smoother_bld.f90 @@ -71,13 +71,12 @@ subroutine mld_z_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold) nrow_a = a%get_nrows() nztota = a%get_nzeros() - if( sm%checkres ) sm%pa => a + sm%pa => a select type (smsv => sm%sv) class is (mld_z_diag_solver_type) call sm%nd%free() - sm%pa => a - sm%nd_nnz_tot = nztota + sm%nd_nnz_tot = nztota call psb_sum(ictxt,sm%nd_nnz_tot) call sm%sv%build(a,desc_a,info,amold=amold,vmold=vmold) diff --git a/mlprec/impl/smoother/mld_z_jac_smoother_clone.f90 b/mlprec/impl/smoother/mld_z_jac_smoother_clone.f90 index 5e2b54f8..3347b05c 100644 --- a/mlprec/impl/smoother/mld_z_jac_smoother_clone.f90 +++ b/mlprec/impl/smoother/mld_z_jac_smoother_clone.f90 @@ -77,6 +77,7 @@ subroutine mld_z_jac_smoother_clone(sm,smout,info) allocate(smout%sv,mold=sm%sv,stat=info) if (info == psb_success_) call sm%sv%clone(smo%sv,info) end if + smo%pa => sm%pa class default info = psb_err_internal_error_ diff --git a/mlprec/mld_c_as_smoother.f90 b/mlprec/mld_c_as_smoother.f90 index 017e1520..c7240bb6 100644 --- a/mlprec/mld_c_as_smoother.f90 +++ b/mlprec/mld_c_as_smoother.f90 @@ -65,6 +65,7 @@ module mld_c_as_smoother ! parent type. ! class(mld_c_base_solver_type), allocatable :: sv ! + type(psb_cspmat_type), pointer :: pa => null() type(psb_cspmat_type) :: nd type(psb_desc_type) :: desc_data integer(psb_ipk_) :: novr, restr, prol diff --git a/mlprec/mld_d_as_smoother.f90 b/mlprec/mld_d_as_smoother.f90 index a560706a..0d942799 100644 --- a/mlprec/mld_d_as_smoother.f90 +++ b/mlprec/mld_d_as_smoother.f90 @@ -65,6 +65,7 @@ module mld_d_as_smoother ! parent type. ! class(mld_d_base_solver_type), allocatable :: sv ! + type(psb_dspmat_type), pointer :: pa => null() type(psb_dspmat_type) :: nd type(psb_desc_type) :: desc_data integer(psb_ipk_) :: novr, restr, prol diff --git a/mlprec/mld_s_as_smoother.f90 b/mlprec/mld_s_as_smoother.f90 index 318cb72d..afaab745 100644 --- a/mlprec/mld_s_as_smoother.f90 +++ b/mlprec/mld_s_as_smoother.f90 @@ -65,6 +65,7 @@ module mld_s_as_smoother ! parent type. ! class(mld_s_base_solver_type), allocatable :: sv ! + type(psb_sspmat_type), pointer :: pa => null() type(psb_sspmat_type) :: nd type(psb_desc_type) :: desc_data integer(psb_ipk_) :: novr, restr, prol diff --git a/mlprec/mld_z_as_smoother.f90 b/mlprec/mld_z_as_smoother.f90 index 146dcb9e..ac8c6517 100644 --- a/mlprec/mld_z_as_smoother.f90 +++ b/mlprec/mld_z_as_smoother.f90 @@ -65,6 +65,7 @@ module mld_z_as_smoother ! parent type. ! class(mld_z_base_solver_type), allocatable :: sv ! + type(psb_zspmat_type), pointer :: pa => null() type(psb_zspmat_type) :: nd type(psb_desc_type) :: desc_data integer(psb_ipk_) :: novr, restr, prol