diff --git a/mlprec/impl/smoother/mld_c_as_smoother_apply.f90 b/mlprec/impl/smoother/mld_c_as_smoother_apply.f90 index b2687cce..a48a90b9 100644 --- a/mlprec/impl/smoother/mld_c_as_smoother_apply.f90 +++ b/mlprec/impl/smoother/mld_c_as_smoother_apply.f90 @@ -105,116 +105,218 @@ subroutine mld_c_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,& end if endif - if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then - ! - ! Shortcut: in this case there is nothing else to be done. - ! - 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 + if (.false.) then + 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) + + ! + ! 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) + ! + select case (init_) + case('Z') + call psb_geaxpby(cone,x,czero,tx,desc_data,info) + if (info == 0) call sm%apply_restr(tx,trans_,aux,info) + call sm%sv%apply(cone,tx,czero,ty,sm%desc_data,trans_,aux,info,init='Z') + + case('Y') + call psb_geaxpby(cone,x,czero,tx,desc_data,info) + if (info == 0) call psb_spmm(-cone,sm%pa,y,cone,ww,desc_data,info,& + & work=aux,trans=trans_) + if (info == 0) call sm%apply_restr(tx,trans_,aux,info) + call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Y') - 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) - - ! - ! 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 sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Z') - - case('Y') - 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 + 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,ww,desc_data,info) + call psb_geaxpby(cone,initu,czero,ty,desc_data,info) + if (info == 0) call psb_spmm(-cone,sm%pa,ty,cone,ww,desc_data,info,& + & work=aux,trans=trans_) + if (info == 0) call sm%apply_restr(ww,trans_,aux,info) + 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 (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='missing initu to smoother_apply') + & a_err='Error in sub_aply Jacobi Sweeps = 1') goto 9999 + endif + + do i=1, sweeps-1 + ! + ! Compute Y(j+1) = Y(j)+D^(-1)*(X-A*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. + ! + if (info == 0) call psb_geaxpby(cone,x,czero,ww,desc_data,info) + if (info == 0) call psb_spmm(-cone,sm%pa,ty,cone,ww,desc_data,info,& + & work=aux,trans=trans_) + + if (info /= psb_success_) exit + if (info == 0) call sm%apply_restr(ww,trans_,aux,info) + call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y') + + if (info /= psb_success_) exit + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + + 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 - 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 (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') + + ! + ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) + ! + call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + + 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_ + call psb_errpush(info,name,& + & i_err=(/itwo,sweeps,izero,izero,izero/)) goto 9999 + endif - - do i=1, sweeps-1 + else + if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then + ! + ! Shortcut: in this case there is nothing else to be done. + ! + 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 + ! + ! + ! Apply multiple sweeps of an AS solver + ! to compute an approximate solution of a linear system. ! - ! 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_geasb(tx,sm%desc_data,info) + call psb_geasb(ty,sm%desc_data,info) + call psb_geasb(ww,sm%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) - if (info == 0) 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 (init_) + 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,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,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 (info == 0) call sm%apply_prol(ty,trans_,aux,info) - - end do - if (info /= psb_success_) then - info=psb_err_internal_error_ + 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. + ! + 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_) + + 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 + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + + 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 + + ! + ! 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,& - & 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 + & i_err=(/itwo,sweeps,izero,izero,izero/)) + goto 9999 - endif + endif + end if - if (.not.(4*isz <= size(work))) then deallocate(aux,stat=info) endif 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 55b0d94a..2cee1718 100644 --- a/mlprec/impl/smoother/mld_c_as_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_c_as_smoother_apply_vect.f90 @@ -107,89 +107,140 @@ subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& end if endif - if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then - ! - ! Shortcut: in this case there is nothing else to be done. - ! - call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,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 - ! - ! - ! Apply multiple sweeps of an AS solver - ! to compute an approximate solution of a linear system. - ! - ! - if (size(wv) < 3) then + if (.true.) then + + if (sweeps > 0) then + ! + ! + ! Apply multiple sweeps of an AS solver + ! to compute an approximate solution of a linear system. + ! + ! + if (size(wv) < 3) then call psb_errpush(psb_err_internal_error_,name,& & a_err='invalid wv size in smoother_apply') - goto 9999 - end if - - ! - ! This is tricky. This smoother has a descriptor sm%desc_data - ! for an index space potentially different from - ! that of desc_data. Hence the size of the work vectors - ! could be wrong. We need to check and reallocate as needed. - ! - do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.& - & (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.& - & (wv(3)%get_nrows() < sm%desc_data%get_local_cols()) - - if (do_realloc_wv) then - call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v) - call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v) - call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v) - end if - - associate(tx => wv(1), ty => wv(2), ww => wv(3)) - - ! Need to zero tx because of the apply_restr call. - call tx%zero() + goto 9999 + end if + ! - ! 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) + ! This is tricky. This smoother has a descriptor sm%desc_data + ! for an index space potentially different from + ! that of desc_data. Hence the size of the work vectors + ! could be wrong. We need to check and reallocate as needed. ! - 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 sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z') - - case('Y') - 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,wv(4:),info,init='Y') - - case('U') - if (.not.present(initu)) then + do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.& + & (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.& + & (wv(3)%get_nrows() < sm%desc_data%get_local_cols()) + + if (do_realloc_wv) then + call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v) + call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v) + call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v) + end if + + associate(tx => wv(1), ty => wv(2), ww => wv(3)) + + ! Need to zero tx because of the apply_restr call. + call tx%zero() + ! + ! 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) + ! + + select case (init_) + case('Z') + call psb_geaxpby(cone,x,czero,tx,desc_data,info) + if (info == 0) call sm%apply_restr(tx,trans_,aux,info) + call sm%sv%apply(cone,tx,czero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z') + + case('Y') + call psb_geaxpby(cone,x,czero,ww,desc_data,info) + if (info == 0) call psb_spmm(-cone,sm%pa,y,cone,ww,desc_data,info,& + & work=aux,trans=trans_) + if (info == 0) call sm%apply_restr(ww,trans_,aux,info) + call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,wv(4:),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,ww,desc_data,info) + call psb_geaxpby(cone,initu,czero,ty,desc_data,info) + if (info == 0) call psb_spmm(-cone,sm%pa,ty,cone,ww,desc_data,info,& + & work=aux,trans=trans_) + if (info == 0) call sm%apply_restr(ww,trans_,aux,info) + call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,wv(4:),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 (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='missing initu to smoother_apply') + & a_err='Error in sub_aply Jacobi Sweeps = 1') goto 9999 + endif + + do i=1, sweeps-1 + ! + ! Compute Y(j+1) = Y(j)+D^(-1)*(X-A*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. + ! + if (info == 0) call psb_geaxpby(cone,x,czero,ww,desc_data,info) + if (info == 0) call psb_spmm(-cone,sm%pa,ty,cone,ww,desc_data,info,& + & work=aux,trans=trans_) + + if (info /= psb_success_) exit + if (info == 0) call sm%apply_restr(ww,trans_,aux,info) + call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y') + + if (info /= psb_success_) exit + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + + 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 - 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,wv(4:),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 (info == 0) call sm%apply_prol(ty,trans_,aux,info) + ! + ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) + ! + call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + end associate + + 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_ + call psb_errpush(info,name,& + & i_err=(/itwo,sweeps,izero,izero,izero/)) + goto 9999 + + endif + + else + + if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then + ! + ! Shortcut: in this case there is nothing else to be done. + ! + call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -197,47 +248,125 @@ subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& goto 9999 endif - do i=1, sweeps-1 + else if (sweeps >= 0) then + ! + ! + ! Apply multiple sweeps of an AS solver + ! to compute an approximate solution of a linear system. + ! + ! + if (size(wv) < 3) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='invalid wv size in smoother_apply') + goto 9999 + end if + + ! + ! This is tricky. This smoother has a descriptor sm%desc_data + ! for an index space potentially different from + ! that of desc_data. Hence the size of the work vectors + ! could be wrong. We need to check and reallocate as needed. + ! + do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.& + & (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.& + & (wv(3)%get_nrows() < sm%desc_data%get_local_cols()) + + if (do_realloc_wv) then + call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v) + call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v) + call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v) + end if + + associate(tx => wv(1), ty => wv(2), ww => wv(3)) + + ! Need to zero tx because of the apply_restr call. + call tx%zero() ! - ! 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. + ! 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) - if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,& - & work=aux,trans=trans_) - if (info /= psb_success_) exit + select case (init_) + case('Z') + call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z') - call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y') + case('Y') + 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,wv(4:),info,init='Y') - if (info /= psb_success_) exit + 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,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,wv(4:),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 (info == 0) call sm%apply_prol(ty,trans_,aux,info) - end do + 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 - 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 + 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. + ! + 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_) - ! - ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) - ! - call psb_geaxpby(alpha,ty,beta,y,desc_data,info) - end associate + if (info /= psb_success_) exit - else - - info = psb_err_iarg_neg_ - call psb_errpush(info,name,& - & i_err=(/itwo,sweeps,izero,izero,izero/)) - goto 9999 + call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y') - endif + if (info /= psb_success_) exit + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + 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 + + ! + ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) + ! + call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + end associate + + else + + info = psb_err_iarg_neg_ + call psb_errpush(info,name,& + & i_err=(/itwo,sweeps,izero,izero,izero/)) + goto 9999 + + endif + end if if (.not.(4*isz <= size(work))) then deallocate(aux,stat=info) diff --git a/mlprec/impl/smoother/mld_d_as_smoother_apply.f90 b/mlprec/impl/smoother/mld_d_as_smoother_apply.f90 index 4a8835d9..11187722 100644 --- a/mlprec/impl/smoother/mld_d_as_smoother_apply.f90 +++ b/mlprec/impl/smoother/mld_d_as_smoother_apply.f90 @@ -105,116 +105,218 @@ subroutine mld_d_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,& end if endif - if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then - ! - ! Shortcut: in this case there is nothing else to be done. - ! - 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 + if (.false.) then + 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) + + ! + ! 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) + ! + select case (init_) + case('Z') + call psb_geaxpby(done,x,dzero,tx,desc_data,info) + if (info == 0) call sm%apply_restr(tx,trans_,aux,info) + call sm%sv%apply(done,tx,dzero,ty,sm%desc_data,trans_,aux,info,init='Z') + + case('Y') + call psb_geaxpby(done,x,dzero,tx,desc_data,info) + if (info == 0) call psb_spmm(-done,sm%pa,y,done,ww,desc_data,info,& + & work=aux,trans=trans_) + if (info == 0) call sm%apply_restr(tx,trans_,aux,info) + call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Y') - 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) - - ! - ! 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 sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Z') - - case('Y') - 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 + 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,ww,desc_data,info) + call psb_geaxpby(done,initu,dzero,ty,desc_data,info) + if (info == 0) call psb_spmm(-done,sm%pa,ty,done,ww,desc_data,info,& + & work=aux,trans=trans_) + if (info == 0) call sm%apply_restr(ww,trans_,aux,info) + 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 (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='missing initu to smoother_apply') + & a_err='Error in sub_aply Jacobi Sweeps = 1') goto 9999 + endif + + do i=1, sweeps-1 + ! + ! Compute Y(j+1) = Y(j)+D^(-1)*(X-A*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. + ! + if (info == 0) call psb_geaxpby(done,x,dzero,ww,desc_data,info) + if (info == 0) call psb_spmm(-done,sm%pa,ty,done,ww,desc_data,info,& + & work=aux,trans=trans_) + + if (info /= psb_success_) exit + if (info == 0) call sm%apply_restr(ww,trans_,aux,info) + call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y') + + if (info /= psb_success_) exit + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + + 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 - 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 (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') + + ! + ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) + ! + call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + + 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_ + call psb_errpush(info,name,& + & i_err=(/itwo,sweeps,izero,izero,izero/)) goto 9999 + endif - - do i=1, sweeps-1 + else + if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then + ! + ! Shortcut: in this case there is nothing else to be done. + ! + 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 + ! + ! + ! Apply multiple sweeps of an AS solver + ! to compute an approximate solution of a linear system. ! - ! 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_geasb(tx,sm%desc_data,info) + call psb_geasb(ty,sm%desc_data,info) + call psb_geasb(ww,sm%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) - if (info == 0) 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 (init_) + 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,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,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 (info == 0) call sm%apply_prol(ty,trans_,aux,info) - - end do - if (info /= psb_success_) then - info=psb_err_internal_error_ + 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. + ! + 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_) + + 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 + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + + 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 + + ! + ! 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,& - & 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 + & i_err=(/itwo,sweeps,izero,izero,izero/)) + goto 9999 - endif + endif + end if - if (.not.(4*isz <= size(work))) then deallocate(aux,stat=info) endif 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 7ae8e8bf..c6f7ce2d 100644 --- a/mlprec/impl/smoother/mld_d_as_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_d_as_smoother_apply_vect.f90 @@ -107,89 +107,140 @@ subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& end if endif - if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then - ! - ! Shortcut: in this case there is nothing else to be done. - ! - call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,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 - ! - ! - ! Apply multiple sweeps of an AS solver - ! to compute an approximate solution of a linear system. - ! - ! - if (size(wv) < 3) then + if (.true.) then + + if (sweeps > 0) then + ! + ! + ! Apply multiple sweeps of an AS solver + ! to compute an approximate solution of a linear system. + ! + ! + if (size(wv) < 3) then call psb_errpush(psb_err_internal_error_,name,& & a_err='invalid wv size in smoother_apply') - goto 9999 - end if - - ! - ! This is tricky. This smoother has a descriptor sm%desc_data - ! for an index space potentially different from - ! that of desc_data. Hence the size of the work vectors - ! could be wrong. We need to check and reallocate as needed. - ! - do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.& - & (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.& - & (wv(3)%get_nrows() < sm%desc_data%get_local_cols()) - - if (do_realloc_wv) then - call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v) - call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v) - call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v) - end if - - associate(tx => wv(1), ty => wv(2), ww => wv(3)) - - ! Need to zero tx because of the apply_restr call. - call tx%zero() + goto 9999 + end if + ! - ! 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) + ! This is tricky. This smoother has a descriptor sm%desc_data + ! for an index space potentially different from + ! that of desc_data. Hence the size of the work vectors + ! could be wrong. We need to check and reallocate as needed. ! - 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 sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z') - - case('Y') - 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,wv(4:),info,init='Y') - - case('U') - if (.not.present(initu)) then + do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.& + & (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.& + & (wv(3)%get_nrows() < sm%desc_data%get_local_cols()) + + if (do_realloc_wv) then + call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v) + call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v) + call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v) + end if + + associate(tx => wv(1), ty => wv(2), ww => wv(3)) + + ! Need to zero tx because of the apply_restr call. + call tx%zero() + ! + ! 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) + ! + + select case (init_) + case('Z') + call psb_geaxpby(done,x,dzero,tx,desc_data,info) + if (info == 0) call sm%apply_restr(tx,trans_,aux,info) + call sm%sv%apply(done,tx,dzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z') + + case('Y') + call psb_geaxpby(done,x,dzero,ww,desc_data,info) + if (info == 0) call psb_spmm(-done,sm%pa,y,done,ww,desc_data,info,& + & work=aux,trans=trans_) + if (info == 0) call sm%apply_restr(ww,trans_,aux,info) + call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,wv(4:),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,ww,desc_data,info) + call psb_geaxpby(done,initu,dzero,ty,desc_data,info) + if (info == 0) call psb_spmm(-done,sm%pa,ty,done,ww,desc_data,info,& + & work=aux,trans=trans_) + if (info == 0) call sm%apply_restr(ww,trans_,aux,info) + call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,wv(4:),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 (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='missing initu to smoother_apply') + & a_err='Error in sub_aply Jacobi Sweeps = 1') goto 9999 + endif + + do i=1, sweeps-1 + ! + ! Compute Y(j+1) = Y(j)+D^(-1)*(X-A*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. + ! + if (info == 0) call psb_geaxpby(done,x,dzero,ww,desc_data,info) + if (info == 0) call psb_spmm(-done,sm%pa,ty,done,ww,desc_data,info,& + & work=aux,trans=trans_) + + if (info /= psb_success_) exit + if (info == 0) call sm%apply_restr(ww,trans_,aux,info) + call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y') + + if (info /= psb_success_) exit + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + + 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 - 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,wv(4:),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 (info == 0) call sm%apply_prol(ty,trans_,aux,info) + ! + ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) + ! + call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + end associate + + 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_ + call psb_errpush(info,name,& + & i_err=(/itwo,sweeps,izero,izero,izero/)) + goto 9999 + + endif + + else + + if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then + ! + ! Shortcut: in this case there is nothing else to be done. + ! + call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -197,47 +248,125 @@ subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& goto 9999 endif - do i=1, sweeps-1 + else if (sweeps >= 0) then + ! + ! + ! Apply multiple sweeps of an AS solver + ! to compute an approximate solution of a linear system. + ! + ! + if (size(wv) < 3) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='invalid wv size in smoother_apply') + goto 9999 + end if + + ! + ! This is tricky. This smoother has a descriptor sm%desc_data + ! for an index space potentially different from + ! that of desc_data. Hence the size of the work vectors + ! could be wrong. We need to check and reallocate as needed. + ! + do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.& + & (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.& + & (wv(3)%get_nrows() < sm%desc_data%get_local_cols()) + + if (do_realloc_wv) then + call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v) + call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v) + call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v) + end if + + associate(tx => wv(1), ty => wv(2), ww => wv(3)) + + ! Need to zero tx because of the apply_restr call. + call tx%zero() ! - ! 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. + ! 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) - if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,& - & work=aux,trans=trans_) - if (info /= psb_success_) exit + select case (init_) + case('Z') + call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z') - call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y') + case('Y') + 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,wv(4:),info,init='Y') - if (info /= psb_success_) exit + 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,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,wv(4:),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 (info == 0) call sm%apply_prol(ty,trans_,aux,info) - end do + 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 - 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 + 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. + ! + 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_) - ! - ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) - ! - call psb_geaxpby(alpha,ty,beta,y,desc_data,info) - end associate + if (info /= psb_success_) exit - else - - info = psb_err_iarg_neg_ - call psb_errpush(info,name,& - & i_err=(/itwo,sweeps,izero,izero,izero/)) - goto 9999 + call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y') - endif + if (info /= psb_success_) exit + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + 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 + + ! + ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) + ! + call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + end associate + + else + + info = psb_err_iarg_neg_ + call psb_errpush(info,name,& + & i_err=(/itwo,sweeps,izero,izero,izero/)) + goto 9999 + + endif + end if if (.not.(4*isz <= size(work))) then deallocate(aux,stat=info) diff --git a/mlprec/impl/smoother/mld_s_as_smoother_apply.f90 b/mlprec/impl/smoother/mld_s_as_smoother_apply.f90 index 0dbb467e..a8d7a2e2 100644 --- a/mlprec/impl/smoother/mld_s_as_smoother_apply.f90 +++ b/mlprec/impl/smoother/mld_s_as_smoother_apply.f90 @@ -105,116 +105,218 @@ subroutine mld_s_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,& end if endif - if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then - ! - ! Shortcut: in this case there is nothing else to be done. - ! - 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 + if (.false.) then + 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) + + ! + ! 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) + ! + select case (init_) + case('Z') + call psb_geaxpby(sone,x,szero,tx,desc_data,info) + if (info == 0) call sm%apply_restr(tx,trans_,aux,info) + call sm%sv%apply(sone,tx,szero,ty,sm%desc_data,trans_,aux,info,init='Z') + + case('Y') + call psb_geaxpby(sone,x,szero,tx,desc_data,info) + if (info == 0) call psb_spmm(-sone,sm%pa,y,sone,ww,desc_data,info,& + & work=aux,trans=trans_) + if (info == 0) call sm%apply_restr(tx,trans_,aux,info) + call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Y') - 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) - - ! - ! 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 sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Z') - - case('Y') - 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 + 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,ww,desc_data,info) + call psb_geaxpby(sone,initu,szero,ty,desc_data,info) + if (info == 0) call psb_spmm(-sone,sm%pa,ty,sone,ww,desc_data,info,& + & work=aux,trans=trans_) + if (info == 0) call sm%apply_restr(ww,trans_,aux,info) + 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 (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='missing initu to smoother_apply') + & a_err='Error in sub_aply Jacobi Sweeps = 1') goto 9999 + endif + + do i=1, sweeps-1 + ! + ! Compute Y(j+1) = Y(j)+D^(-1)*(X-A*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. + ! + if (info == 0) call psb_geaxpby(sone,x,szero,ww,desc_data,info) + if (info == 0) call psb_spmm(-sone,sm%pa,ty,sone,ww,desc_data,info,& + & work=aux,trans=trans_) + + if (info /= psb_success_) exit + if (info == 0) call sm%apply_restr(ww,trans_,aux,info) + call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y') + + if (info /= psb_success_) exit + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + + 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 - 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 (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') + + ! + ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) + ! + call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + + 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_ + call psb_errpush(info,name,& + & i_err=(/itwo,sweeps,izero,izero,izero/)) goto 9999 + endif - - do i=1, sweeps-1 + else + if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then + ! + ! Shortcut: in this case there is nothing else to be done. + ! + 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 + ! + ! + ! Apply multiple sweeps of an AS solver + ! to compute an approximate solution of a linear system. ! - ! 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_geasb(tx,sm%desc_data,info) + call psb_geasb(ty,sm%desc_data,info) + call psb_geasb(ww,sm%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) - if (info == 0) 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 (init_) + 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,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,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 (info == 0) call sm%apply_prol(ty,trans_,aux,info) - - end do - if (info /= psb_success_) then - info=psb_err_internal_error_ + 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. + ! + 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_) + + 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 + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + + 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 + + ! + ! 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,& - & 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 + & i_err=(/itwo,sweeps,izero,izero,izero/)) + goto 9999 - endif + endif + end if - if (.not.(4*isz <= size(work))) then deallocate(aux,stat=info) endif 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 e20643ba..5ad8c380 100644 --- a/mlprec/impl/smoother/mld_s_as_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_s_as_smoother_apply_vect.f90 @@ -107,89 +107,140 @@ subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& end if endif - if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then - ! - ! Shortcut: in this case there is nothing else to be done. - ! - call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,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 - ! - ! - ! Apply multiple sweeps of an AS solver - ! to compute an approximate solution of a linear system. - ! - ! - if (size(wv) < 3) then + if (.true.) then + + if (sweeps > 0) then + ! + ! + ! Apply multiple sweeps of an AS solver + ! to compute an approximate solution of a linear system. + ! + ! + if (size(wv) < 3) then call psb_errpush(psb_err_internal_error_,name,& & a_err='invalid wv size in smoother_apply') - goto 9999 - end if - - ! - ! This is tricky. This smoother has a descriptor sm%desc_data - ! for an index space potentially different from - ! that of desc_data. Hence the size of the work vectors - ! could be wrong. We need to check and reallocate as needed. - ! - do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.& - & (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.& - & (wv(3)%get_nrows() < sm%desc_data%get_local_cols()) - - if (do_realloc_wv) then - call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v) - call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v) - call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v) - end if - - associate(tx => wv(1), ty => wv(2), ww => wv(3)) - - ! Need to zero tx because of the apply_restr call. - call tx%zero() + goto 9999 + end if + ! - ! 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) + ! This is tricky. This smoother has a descriptor sm%desc_data + ! for an index space potentially different from + ! that of desc_data. Hence the size of the work vectors + ! could be wrong. We need to check and reallocate as needed. ! - 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 sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z') - - case('Y') - 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,wv(4:),info,init='Y') - - case('U') - if (.not.present(initu)) then + do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.& + & (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.& + & (wv(3)%get_nrows() < sm%desc_data%get_local_cols()) + + if (do_realloc_wv) then + call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v) + call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v) + call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v) + end if + + associate(tx => wv(1), ty => wv(2), ww => wv(3)) + + ! Need to zero tx because of the apply_restr call. + call tx%zero() + ! + ! 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) + ! + + select case (init_) + case('Z') + call psb_geaxpby(sone,x,szero,tx,desc_data,info) + if (info == 0) call sm%apply_restr(tx,trans_,aux,info) + call sm%sv%apply(sone,tx,szero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z') + + case('Y') + call psb_geaxpby(sone,x,szero,ww,desc_data,info) + if (info == 0) call psb_spmm(-sone,sm%pa,y,sone,ww,desc_data,info,& + & work=aux,trans=trans_) + if (info == 0) call sm%apply_restr(ww,trans_,aux,info) + call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,wv(4:),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,ww,desc_data,info) + call psb_geaxpby(sone,initu,szero,ty,desc_data,info) + if (info == 0) call psb_spmm(-sone,sm%pa,ty,sone,ww,desc_data,info,& + & work=aux,trans=trans_) + if (info == 0) call sm%apply_restr(ww,trans_,aux,info) + call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,wv(4:),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 (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='missing initu to smoother_apply') + & a_err='Error in sub_aply Jacobi Sweeps = 1') goto 9999 + endif + + do i=1, sweeps-1 + ! + ! Compute Y(j+1) = Y(j)+D^(-1)*(X-A*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. + ! + if (info == 0) call psb_geaxpby(sone,x,szero,ww,desc_data,info) + if (info == 0) call psb_spmm(-sone,sm%pa,ty,sone,ww,desc_data,info,& + & work=aux,trans=trans_) + + if (info /= psb_success_) exit + if (info == 0) call sm%apply_restr(ww,trans_,aux,info) + call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y') + + if (info /= psb_success_) exit + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + + 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 - 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,wv(4:),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 (info == 0) call sm%apply_prol(ty,trans_,aux,info) + ! + ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) + ! + call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + end associate + + 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_ + call psb_errpush(info,name,& + & i_err=(/itwo,sweeps,izero,izero,izero/)) + goto 9999 + + endif + + else + + if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then + ! + ! Shortcut: in this case there is nothing else to be done. + ! + call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -197,47 +248,125 @@ subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& goto 9999 endif - do i=1, sweeps-1 + else if (sweeps >= 0) then + ! + ! + ! Apply multiple sweeps of an AS solver + ! to compute an approximate solution of a linear system. + ! + ! + if (size(wv) < 3) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='invalid wv size in smoother_apply') + goto 9999 + end if + + ! + ! This is tricky. This smoother has a descriptor sm%desc_data + ! for an index space potentially different from + ! that of desc_data. Hence the size of the work vectors + ! could be wrong. We need to check and reallocate as needed. + ! + do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.& + & (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.& + & (wv(3)%get_nrows() < sm%desc_data%get_local_cols()) + + if (do_realloc_wv) then + call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v) + call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v) + call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v) + end if + + associate(tx => wv(1), ty => wv(2), ww => wv(3)) + + ! Need to zero tx because of the apply_restr call. + call tx%zero() ! - ! 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. + ! 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) - if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,& - & work=aux,trans=trans_) - if (info /= psb_success_) exit + select case (init_) + case('Z') + call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z') - call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y') + case('Y') + 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,wv(4:),info,init='Y') - if (info /= psb_success_) exit + 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,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,wv(4:),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 (info == 0) call sm%apply_prol(ty,trans_,aux,info) - end do + 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 - 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 + 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. + ! + 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_) - ! - ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) - ! - call psb_geaxpby(alpha,ty,beta,y,desc_data,info) - end associate + if (info /= psb_success_) exit - else - - info = psb_err_iarg_neg_ - call psb_errpush(info,name,& - & i_err=(/itwo,sweeps,izero,izero,izero/)) - goto 9999 + call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y') - endif + if (info /= psb_success_) exit + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + 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 + + ! + ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) + ! + call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + end associate + + else + + info = psb_err_iarg_neg_ + call psb_errpush(info,name,& + & i_err=(/itwo,sweeps,izero,izero,izero/)) + goto 9999 + + endif + end if if (.not.(4*isz <= size(work))) then deallocate(aux,stat=info) diff --git a/mlprec/impl/smoother/mld_z_as_smoother_apply.f90 b/mlprec/impl/smoother/mld_z_as_smoother_apply.f90 index a0e18c1a..6179bd50 100644 --- a/mlprec/impl/smoother/mld_z_as_smoother_apply.f90 +++ b/mlprec/impl/smoother/mld_z_as_smoother_apply.f90 @@ -105,116 +105,218 @@ subroutine mld_z_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,& end if endif - if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then - ! - ! Shortcut: in this case there is nothing else to be done. - ! - 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 + if (.false.) then + 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) + + ! + ! 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) + ! + select case (init_) + case('Z') + call psb_geaxpby(zone,x,zzero,tx,desc_data,info) + if (info == 0) call sm%apply_restr(tx,trans_,aux,info) + call sm%sv%apply(zone,tx,zzero,ty,sm%desc_data,trans_,aux,info,init='Z') + + case('Y') + call psb_geaxpby(zone,x,zzero,tx,desc_data,info) + if (info == 0) call psb_spmm(-zone,sm%pa,y,zone,ww,desc_data,info,& + & work=aux,trans=trans_) + if (info == 0) call sm%apply_restr(tx,trans_,aux,info) + call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Y') - 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) - - ! - ! 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 sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Z') - - case('Y') - 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 + 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,ww,desc_data,info) + call psb_geaxpby(zone,initu,zzero,ty,desc_data,info) + if (info == 0) call psb_spmm(-zone,sm%pa,ty,zone,ww,desc_data,info,& + & work=aux,trans=trans_) + if (info == 0) call sm%apply_restr(ww,trans_,aux,info) + 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 (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='missing initu to smoother_apply') + & a_err='Error in sub_aply Jacobi Sweeps = 1') goto 9999 + endif + + do i=1, sweeps-1 + ! + ! Compute Y(j+1) = Y(j)+D^(-1)*(X-A*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. + ! + if (info == 0) call psb_geaxpby(zone,x,zzero,ww,desc_data,info) + if (info == 0) call psb_spmm(-zone,sm%pa,ty,zone,ww,desc_data,info,& + & work=aux,trans=trans_) + + if (info /= psb_success_) exit + if (info == 0) call sm%apply_restr(ww,trans_,aux,info) + call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y') + + if (info /= psb_success_) exit + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + + 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 - 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 (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') + + ! + ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) + ! + call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + + 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_ + call psb_errpush(info,name,& + & i_err=(/itwo,sweeps,izero,izero,izero/)) goto 9999 + endif - - do i=1, sweeps-1 + else + if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then + ! + ! Shortcut: in this case there is nothing else to be done. + ! + 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 + ! + ! + ! Apply multiple sweeps of an AS solver + ! to compute an approximate solution of a linear system. ! - ! 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_geasb(tx,sm%desc_data,info) + call psb_geasb(ty,sm%desc_data,info) + call psb_geasb(ww,sm%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) - if (info == 0) 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 (init_) + 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,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,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 (info == 0) call sm%apply_prol(ty,trans_,aux,info) - - end do - if (info /= psb_success_) then - info=psb_err_internal_error_ + 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. + ! + 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_) + + 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 + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + + 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 + + ! + ! 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,& - & 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 + & i_err=(/itwo,sweeps,izero,izero,izero/)) + goto 9999 - endif + endif + end if - if (.not.(4*isz <= size(work))) then deallocate(aux,stat=info) endif 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 53bf90f6..f80fc67b 100644 --- a/mlprec/impl/smoother/mld_z_as_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_z_as_smoother_apply_vect.f90 @@ -107,89 +107,140 @@ subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& end if endif - if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then - ! - ! Shortcut: in this case there is nothing else to be done. - ! - call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,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 - ! - ! - ! Apply multiple sweeps of an AS solver - ! to compute an approximate solution of a linear system. - ! - ! - if (size(wv) < 3) then + if (.true.) then + + if (sweeps > 0) then + ! + ! + ! Apply multiple sweeps of an AS solver + ! to compute an approximate solution of a linear system. + ! + ! + if (size(wv) < 3) then call psb_errpush(psb_err_internal_error_,name,& & a_err='invalid wv size in smoother_apply') - goto 9999 - end if - - ! - ! This is tricky. This smoother has a descriptor sm%desc_data - ! for an index space potentially different from - ! that of desc_data. Hence the size of the work vectors - ! could be wrong. We need to check and reallocate as needed. - ! - do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.& - & (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.& - & (wv(3)%get_nrows() < sm%desc_data%get_local_cols()) - - if (do_realloc_wv) then - call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v) - call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v) - call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v) - end if - - associate(tx => wv(1), ty => wv(2), ww => wv(3)) - - ! Need to zero tx because of the apply_restr call. - call tx%zero() + goto 9999 + end if + ! - ! 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) + ! This is tricky. This smoother has a descriptor sm%desc_data + ! for an index space potentially different from + ! that of desc_data. Hence the size of the work vectors + ! could be wrong. We need to check and reallocate as needed. ! - 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 sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z') - - case('Y') - 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,wv(4:),info,init='Y') - - case('U') - if (.not.present(initu)) then + do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.& + & (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.& + & (wv(3)%get_nrows() < sm%desc_data%get_local_cols()) + + if (do_realloc_wv) then + call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v) + call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v) + call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v) + end if + + associate(tx => wv(1), ty => wv(2), ww => wv(3)) + + ! Need to zero tx because of the apply_restr call. + call tx%zero() + ! + ! 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) + ! + + select case (init_) + case('Z') + call psb_geaxpby(zone,x,zzero,tx,desc_data,info) + if (info == 0) call sm%apply_restr(tx,trans_,aux,info) + call sm%sv%apply(zone,tx,zzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z') + + case('Y') + call psb_geaxpby(zone,x,zzero,ww,desc_data,info) + if (info == 0) call psb_spmm(-zone,sm%pa,y,zone,ww,desc_data,info,& + & work=aux,trans=trans_) + if (info == 0) call sm%apply_restr(ww,trans_,aux,info) + call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,wv(4:),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,ww,desc_data,info) + call psb_geaxpby(zone,initu,zzero,ty,desc_data,info) + if (info == 0) call psb_spmm(-zone,sm%pa,ty,zone,ww,desc_data,info,& + & work=aux,trans=trans_) + if (info == 0) call sm%apply_restr(ww,trans_,aux,info) + call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,wv(4:),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 (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='missing initu to smoother_apply') + & a_err='Error in sub_aply Jacobi Sweeps = 1') goto 9999 + endif + + do i=1, sweeps-1 + ! + ! Compute Y(j+1) = Y(j)+D^(-1)*(X-A*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. + ! + if (info == 0) call psb_geaxpby(zone,x,zzero,ww,desc_data,info) + if (info == 0) call psb_spmm(-zone,sm%pa,ty,zone,ww,desc_data,info,& + & work=aux,trans=trans_) + + if (info /= psb_success_) exit + if (info == 0) call sm%apply_restr(ww,trans_,aux,info) + call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y') + + if (info /= psb_success_) exit + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + + 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 - 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,wv(4:),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 (info == 0) call sm%apply_prol(ty,trans_,aux,info) + ! + ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) + ! + call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + end associate + + 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_ + call psb_errpush(info,name,& + & i_err=(/itwo,sweeps,izero,izero,izero/)) + goto 9999 + + endif + + else + + if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then + ! + ! Shortcut: in this case there is nothing else to be done. + ! + call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -197,47 +248,125 @@ subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& goto 9999 endif - do i=1, sweeps-1 + else if (sweeps >= 0) then + ! + ! + ! Apply multiple sweeps of an AS solver + ! to compute an approximate solution of a linear system. + ! + ! + if (size(wv) < 3) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='invalid wv size in smoother_apply') + goto 9999 + end if + + ! + ! This is tricky. This smoother has a descriptor sm%desc_data + ! for an index space potentially different from + ! that of desc_data. Hence the size of the work vectors + ! could be wrong. We need to check and reallocate as needed. + ! + do_realloc_wv = (wv(1)%get_nrows() < sm%desc_data%get_local_cols()).or.& + & (wv(2)%get_nrows() < sm%desc_data%get_local_cols()).or.& + & (wv(3)%get_nrows() < sm%desc_data%get_local_cols()) + + if (do_realloc_wv) then + call psb_geasb(wv(1),sm%desc_data,info,scratch=.true.,mold=wv(2)%v) + call psb_geasb(wv(2),sm%desc_data,info,scratch=.true.,mold=wv(1)%v) + call psb_geasb(wv(3),sm%desc_data,info,scratch=.true.,mold=wv(1)%v) + end if + + associate(tx => wv(1), ty => wv(2), ww => wv(3)) + + ! Need to zero tx because of the apply_restr call. + call tx%zero() ! - ! 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. + ! 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) - if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,& - & work=aux,trans=trans_) - if (info /= psb_success_) exit + select case (init_) + case('Z') + call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Z') - call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y') + case('Y') + 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,wv(4:),info,init='Y') - if (info /= psb_success_) exit + 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,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,wv(4:),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 (info == 0) call sm%apply_prol(ty,trans_,aux,info) - end do + 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 - 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 + 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. + ! + 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_) - ! - ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) - ! - call psb_geaxpby(alpha,ty,beta,y,desc_data,info) - end associate + if (info /= psb_success_) exit - else - - info = psb_err_iarg_neg_ - call psb_errpush(info,name,& - & i_err=(/itwo,sweeps,izero,izero,izero/)) - goto 9999 + call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,wv(4:),info,init='Y') - endif + if (info /= psb_success_) exit + if (info == 0) call sm%apply_prol(ty,trans_,aux,info) + 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 + + ! + ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) + ! + call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + end associate + + else + + info = psb_err_iarg_neg_ + call psb_errpush(info,name,& + & i_err=(/itwo,sweeps,izero,izero,izero/)) + goto 9999 + + endif + end if if (.not.(4*isz <= size(work))) then deallocate(aux,stat=info)