diff --git a/src/prec/Makefile b/src/prec/Makefile index 28f5d0a8..30ba5736 100644 --- a/src/prec/Makefile +++ b/src/prec/Makefile @@ -6,8 +6,11 @@ LIBDIR=../../lib/ MPFOBJS=psb_dilu_bld.o psb_dbldaggrmat.o F90OBJS=psb_dasmatbld.o psb_dslu_bld.o psb_dumf_bld.o psb_dilu_fct.o\ psb_dmlprc_bld.o psb_dsp_renum.o\ - psb_dprec.o psb_dprecbld.o gps.o psb_dprecfree.o psb_dprecset.o \ - psb_dbaseprc_bld.o psb_dgenaggrmap.o $(MPFOBJS) + psb_dprecbld.o gps.o psb_dprecfree.o psb_dprecset.o \ + psb_dbaseprc_bld.o psb_dgenaggrmap.o \ + psb_dprc_aply.o psb_dmlprc_aply.o \ + psb_dbaseprc_aply.o psb_dbjac_aply.o\ + $(MPFOBJS) COBJS=psb_slu_impl.o psb_umf_impl.o INCDIRS=-I. -I.. -I$(LIBDIR) diff --git a/src/prec/psb_dbaseprc_aply.f90 b/src/prec/psb_dbaseprc_aply.f90 new file mode 100644 index 00000000..5d839a2b --- /dev/null +++ b/src/prec/psb_dbaseprc_aply.f90 @@ -0,0 +1,296 @@ + +!!$ +!!$ +!!$ MPcube: Multilevel Parallel Preconditioners Package +!!$ for +!!$ Parallel Sparse BLAS v2.0 +!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ Daniela Di Serafino II University of Naples +!!$ Pasqua D'Ambra ICAR-CNR +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MPCUBE group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MPCUBE GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ + +subroutine psb_dbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info) + ! + ! Compute Y <- beta*Y + K^-1 X + ! where K is a a basic preconditioner stored in prec + ! + + use psb_serial_mod + use psb_descriptor_type + use psb_prec_type + use psb_psblas_mod + use psb_const_mod + use psb_error_mod + implicit none + + type(psb_desc_type),intent(in) :: desc_data + type(psb_dbaseprc_type), intent(in) :: prec + real(kind(0.d0)),intent(inout) :: x(:), y(:) + real(kind(0.d0)),intent(in) :: beta + character(len=1) :: trans + real(kind(0.d0)),target :: work(:) + integer, intent(out) :: info + + ! Local variables + integer :: n_row,n_col, int_err(5) + real(kind(1.d0)), pointer :: ww(:), aux(:), tx(:),ty(:) + character ::diagl, diagu + integer :: icontxt,nprow,npcol,me,mycol,i, isz, nrg, err_act + real(kind(1.d0)) :: t1, t2, t3, t4, t5, t6, t7, mpi_wtime + logical,parameter :: debug=.false., debugprt=.false. + real(kind(1.d0)), parameter :: one=1.d0, zero=0.d0 + external mpi_wtime + character(len=20) :: name, ch_err + + interface psb_bjac_aply + subroutine psb_dbjac_aply(prec,x,beta,y,desc_data,trans,work,info) + use psb_descriptor_type + use psb_prec_type + type(psb_desc_type), intent(in) :: desc_data + type(psb_dbaseprc_type), intent(in) :: prec + real(kind(0.d0)),intent(inout) :: x(:), y(:) + real(kind(0.d0)),intent(in) :: beta + character(len=1) :: trans + real(kind(0.d0)),target :: work(:) + integer, intent(out) :: info + end subroutine psb_dbjac_aply + end interface + + name='psb_dbaseprc_aply' + info = 0 + call psb_erractionsave(err_act) + + icontxt=desc_data%matrix_data(psb_ctxt_) + call blacs_gridinfo(icontxt,nprow,npcol,me,mycol) + + diagl='U' + diagu='U' + + select case(trans) + case('N','n') + case('T','t','C','c') + case default + info=40 + int_err(1)=6 + ch_err(2:2)=trans + goto 9999 + end select + + select case(prec%iprcparm(p_type_)) + + case(noprec_) + + n_row=desc_data%matrix_data(psb_n_row_) + if (beta==zero) then + y(1:n_row) = x(1:n_row) + else if (beta==one) then + y(1:n_row) = x(1:n_row) + y(1:n_row) + else if (beta==-one) then + y(1:n_row) = x(1:n_row) - y(1:n_row) + else + y(1:n_row) = x(1:n_row) + beta*y(1:n_row) + end if + + case(diagsc_) + + n_row=desc_data%matrix_data(psb_n_row_) + if (beta==zero) then + y(1:n_row) = x(1:n_row)*prec%d(1:n_row) + else if (beta==one) then + y(1:n_row) = x(1:n_row)*prec%d(1:n_row) + y(1:n_row) + else if (beta==-one) then + y(1:n_row) = x(1:n_row)*prec%d(1:n_row) - y(1:n_row) + else + y(1:n_row) = x(1:n_row)*prec%d(1:n_row) + beta*y(1:n_row) + end if + + case(bja_) + + call psb_bjac_aply(prec,x,beta,y,desc_data,trans,work,info) + if(info.ne.0) then + info=4010 + ch_err='psb_bjac_aply' + goto 9999 + end if + + case(asm_,ras_,ash_,rash_) + + if (prec%iprcparm(n_ovr_)==0) then + ! shortcut: this fixes performance for RAS(0) == BJA + call psb_bjac_aply(prec,x,beta,y,desc_data,trans,work,info) + if(info.ne.0) then + info=4010 + ch_err='psb_bjacaply' + goto 9999 + end if + + else + ! Note: currently trans is unused. + n_row=prec%desc_data%matrix_data(psb_n_row_) + n_col=prec%desc_data%matrix_data(psb_n_col_) + + isz=max(n_row,N_COL) + if ((6*isz) <= size(work)) then + ww => work(1:isz) + tx => work(isz+1:2*isz) + ty => work(2*isz+1:3*isz) + aux => work(3*isz+1:) + else if ((4*isz) <= size(work)) then + aux => work(1:) + allocate(ww(isz),tx(isz),ty(isz),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + else if ((3*isz) <= size(work)) then + ww => work(1:isz) + tx => work(isz+1:2*isz) + ty => work(2*isz+1:3*isz) + allocate(aux(4*isz),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + + else + allocate(ww(isz),tx(isz),ty(isz),& + &aux(4*isz),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + + endif + + if (debugprt) write(0,*)' vdiag: ',prec%d(:) + if (debug) write(0,*) 'Bi-CGSTAB with Additive Schwarz prec' + + tx(1:desc_data%matrix_data(psb_n_row_)) = x(1:desc_data%matrix_data(psb_n_row_)) + tx(desc_data%matrix_data(psb_n_row_)+1:isz) = zero + + if (prec%iprcparm(restr_)==psb_halo_) then + call psb_halo(tx,prec%desc_data,info,work=aux) + if(info /=0) then + info=4010 + ch_err='psb_halo' + goto 9999 + end if + else if (prec%iprcparm(restr_) /= psb_none_) then + write(0,*) 'Problem in PRC_APLY: Unknown value for restriction ',& + &prec%iprcparm(restr_) + end if + + if (prec%iprcparm(iren_)>0) then + call dgelp('N',n_row,1,prec%perm,tx,isz,ww,isz,info) + if(info /=0) then + info=4010 + ch_err='psb_dgelp' + goto 9999 + end if + endif + + call psb_bjac_aply(prec,tx,zero,ty,prec%desc_data,trans,aux,info) + if(info.ne.0) then + info=4010 + ch_err='psb_bjac_aply' + goto 9999 + end if + + if (prec%iprcparm(iren_)>0) then + call dgelp('N',n_row,1,prec%invperm,ty,isz,ww,isz,info) + if(info /=0) then + info=4010 + ch_err='psb_dgelp' + goto 9999 + end if + endif + + select case (prec%iprcparm(prol_)) + + case(psb_none_) + ! Would work anyway, but since it's supposed to do nothing... + ! call f90_psovrl(ty,prec%desc_data,update_type=prec%a_restrict) + + case(psb_sum_,psb_avg_) + call psb_ovrl(ty,prec%desc_data,info,& + & update_type=prec%iprcparm(prol_),work=aux) + if(info /=0) then + info=4010 + ch_err='psb_ovrl' + goto 9999 + end if + + case default + write(0,*) 'Problem in PRC_APLY: Unknown value for prolongation ',& + & prec%iprcparm(prol_) + end select + + if (beta == zero) then + y(1:desc_data%matrix_data(psb_n_row_)) = ty(1:desc_data%matrix_data(psb_n_row_)) + else if (beta == one) then + y(1:desc_data%matrix_data(psb_n_row_)) = y(1:desc_data%matrix_data(psb_n_row_)) + & + & ty(1:desc_data%matrix_data(psb_n_row_)) + else if (beta == -one) then + y(1:desc_data%matrix_data(psb_n_row_)) = -y(1:desc_data%matrix_data(psb_n_row_)) + & + & ty(1:desc_data%matrix_data(psb_n_row_)) + else + y(1:desc_data%matrix_data(psb_n_row_)) = beta*y(1:desc_data%matrix_data(psb_n_row_)) + & + & ty(1:desc_data%matrix_data(psb_n_row_)) + end if + + + if ((6*isz) <= size(work)) then + else if ((4*isz) <= size(work)) then + deallocate(ww,tx,ty) + else if ((3*isz) <= size(work)) then + deallocate(aux) + else + deallocate(ww,aux,tx,ty) + endif + end if + case default + write(0,*) 'Invalid PRE%PREC ',prec%iprcparm(p_type_),':',& + & min_prec_,noprec_,diagsc_,bja_,& + & asm_,ras_,ash_,rash_ + end select + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_errpush(info,name,i_err=int_err,a_err=ch_err) + call psb_erractionrestore(err_act) + if (err_act.eq.act_abort) then + call psb_error() + return + end if + return + +end subroutine psb_dbaseprc_aply + diff --git a/src/prec/psb_dbjac_aply.f90 b/src/prec/psb_dbjac_aply.f90 new file mode 100644 index 00000000..70c5790d --- /dev/null +++ b/src/prec/psb_dbjac_aply.f90 @@ -0,0 +1,298 @@ + +!!$ +!!$ +!!$ MPcube: Multilevel Parallel Preconditioners Package +!!$ for +!!$ Parallel Sparse BLAS v2.0 +!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ Daniela Di Serafino II University of Naples +!!$ Pasqua D'Ambra ICAR-CNR +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MPCUBE group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MPCUBE GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +subroutine psb_dbjac_aply(prec,x,beta,y,desc_data,trans,work,info) + ! + ! Compute Y <- beta*Y + K^-1 X + ! where K is a a Block Jacobi preconditioner stored in prec + ! Note that desc_data may or may not be the same as prec%desc_data, + ! but since both are INTENT(IN) this should be legal. + ! + + use psb_serial_mod + use psb_descriptor_type + use psb_prec_type + use psb_psblas_mod + use psb_const_mod + use psb_error_mod + implicit none + + type(psb_desc_type), intent(in) :: desc_data + type(psb_dbaseprc_type), intent(in) :: prec + real(kind(0.d0)),intent(inout) :: x(:), y(:) + real(kind(0.d0)),intent(in) :: beta + character(len=1) :: trans + real(kind(0.d0)),target :: work(:) + integer, intent(out) :: info + + ! Local variables + integer :: n_row,n_col + real(kind(1.d0)), pointer :: ww(:), aux(:), tx(:),ty(:),tb(:) + character ::diagl, diagu + integer :: icontxt,nprow,npcol,me,mycol,i, isz, nrg, err_act, int_err(5) + real(kind(1.d0)) :: t1, t2, t3, t4, t5, t6, t7, mpi_wtime + logical,parameter :: debug=.false., debugprt=.false. + real(kind(1.d0)), parameter :: one=1.d0, zero=0.d0 + external mpi_wtime + character(len=20) :: name, ch_err + + name='psb_bjac_aply' + info = 0 + call psb_erractionsave(err_act) + + icontxt=desc_data%matrix_data(psb_ctxt_) + call blacs_gridinfo(icontxt,nprow,npcol,me,mycol) + + diagl='U' + diagu='U' + + select case(trans) + case('N','n') + case('T','t','C','c') + case default + call psb_errpush(40,name) + goto 9999 + end select + + + n_row=desc_data%matrix_data(psb_n_row_) + n_col=desc_data%matrix_data(psb_n_col_) + + if (n_col <= size(work)) then + ww => work(1:n_col) + if ((4*n_col+n_col) <= size(work)) then + aux => work(n_col+1:) + else + allocate(aux(4*n_col),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + + endif + else + allocate(ww(n_col),aux(4*n_col),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + endif + + + if (prec%iprcparm(jac_sweeps_) == 1) then + + + select case(prec%iprcparm(f_type_)) + case(f_ilu_n_,f_ilu_e_) + + select case(trans) + case('N','n') + + call psb_spsm(one,prec%av(l_pr_),x,zero,ww,desc_data,info,& + & trans='N',unit=diagl,choice=psb_none_,work=aux) + if(info /=0) goto 9999 + ww(1:n_row) = ww(1:n_row)*prec%d(1:n_row) + call psb_spsm(one,prec%av(u_pr_),ww,beta,y,desc_data,info,& + & trans='N',unit=diagu,choice=psb_none_, work=aux) + if(info /=0) goto 9999 + + case('T','t','C','c') + call psb_spsm(one,prec%av(u_pr_),x,zero,ww,desc_data,info,& + & trans=trans,unit=diagu,choice=psb_none_, work=aux) + if(info /=0) goto 9999 + ww(1:n_row) = ww(1:n_row)*prec%d(1:n_row) + call psb_spsm(one,prec%av(l_pr_),ww,beta,y,desc_data,info,& + & trans=trans,unit=diagl,choice=psb_none_,work=aux) + if(info /=0) goto 9999 + + end select + + case(f_slu_) + + ww(1:n_row) = x(1:n_row) + + select case(trans) + case('N','n') + call psb_slu_solve(0,n_row,1,ww,n_row,prec%iprcparm(slu_ptr_),info) + case('T','t','C','c') + call psb_slu_solve(1,n_row,1,ww,n_row,prec%iprcparm(slu_ptr_),info) + end select + + if(info /=0) goto 9999 + + if (beta == 0.d0) then + y(1:n_row) = ww(1:n_row) + else if (beta==1.d0) then + y(1:n_row) = ww(1:n_row) + y(1:n_row) + else if (beta==-1.d0) then + y(1:n_row) = ww(1:n_row) - y(1:n_row) + else + y(1:n_row) = ww(1:n_row) + beta*y(1:n_row) + endif + case (f_umf_) + + + select case(trans) + case('N','n') + call psb_umf_solve(0,n_row,ww,x,n_row,prec%iprcparm(umf_numptr_),info) + case('T','t','C','c') + call psb_umf_solve(1,n_row,ww,x,n_row,prec%iprcparm(umf_numptr_),info) + end select + + if(info /=0) goto 9999 + + if (beta == 0.d0) then + y(1:n_row) = ww(1:n_row) + else if (beta==1.d0) then + y(1:n_row) = ww(1:n_row) + y(1:n_row) + else if (beta==-1.d0) then + y(1:n_row) = ww(1:n_row) - y(1:n_row) + else + y(1:n_row) = ww(1:n_row) + beta*y(1:n_row) + endif + + case default + write(0,*) 'Unknown factorization type in bjac_aply',prec%iprcparm(f_type_) + end select + if (debugprt) write(0,*)' Y: ',y(:) + + else if (prec%iprcparm(jac_sweeps_) > 1) then + + ! Note: we have to add TRANS to this one !!!!!!!!! + + if (size(prec%av) < ap_nd_) then + info = 4011 + goto 9999 + endif + + allocate(tx(n_col),ty(n_col),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + + tx = zero + ty = zero + select case(prec%iprcparm(f_type_)) + case(f_ilu_n_,f_ilu_e_) + do i=1, prec%iprcparm(jac_sweeps_) + ! X(k+1) = M^-1*(b-N*X(k)) + ty(1:n_row) = x(1:n_row) + call psb_spmm(-one,prec%av(ap_nd_),tx,one,ty,& + & prec%desc_data,info,work=aux) + if(info /=0) goto 9999 + call psb_spsm(one,prec%av(l_pr_),ty,zero,ww,& + & prec%desc_data,info,& + & trans='N',unit='U',choice=psb_none_,work=aux) + if(info /=0) goto 9999 + ww(1:n_row) = ww(1:n_row)*prec%d(1:n_row) + call psb_spsm(one,prec%av(u_pr_),ww,zero,tx,& + & prec%desc_data,info,& + & trans='N',unit='U',choice=psb_none_,work=aux) + if(info /=0) goto 9999 + end do + + case(f_slu_) + do i=1, prec%iprcparm(jac_sweeps_) + ! X(k+1) = M^-1*(b-N*X(k)) + ty(1:n_row) = x(1:n_row) + call psb_spmm(-one,prec%av(ap_nd_),tx,one,ty,& + & prec%desc_data,info,work=aux) + if(info /=0) goto 9999 + + call psb_slu_solve(0,n_row,1,ty,n_row,prec%iprcparm(slu_ptr_),info) + if(info /=0) goto 9999 + tx(1:n_row) = ty(1:n_row) + end do + case(f_umf_) + do i=1, prec%iprcparm(jac_sweeps_) + ! X(k+1) = M^-1*(b-N*X(k)) + ty(1:n_row) = x(1:n_row) + call psb_spmm(-one,prec%av(ap_nd_),tx,one,ty,& + & prec%desc_data,info,work=aux) + if(info /=0) goto 9999 + + call psb_umf_solve(0,n_row,ww,ty,n_row,& + & prec%iprcparm(umf_numptr_),info) + if(info /=0) goto 9999 + tx(1:n_row) = ww(1:n_row) + end do + + end select + + if (beta == 0.d0) then + y(1:n_row) = tx(1:n_row) + else if (beta==1.d0) then + y(1:n_row) = tx(1:n_row) + y(1:n_row) + else if (beta==-1.d0) then + y(1:n_row) = tx(1:n_row) - y(1:n_row) + else + y(1:n_row) = tx(1:n_row) + beta*y(1:n_row) + endif + + deallocate(tx,ty) + + + else + + goto 9999 + + endif + + if (n_col <= size(work)) then + if ((4*n_col+n_col) <= size(work)) then + else + deallocate(aux) + endif + else + deallocate(ww,aux) + endif + + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_errpush(info,name,i_err=int_err,a_err=ch_err) + call psb_erractionrestore(err_act) + if (err_act.eq.act_abort) then + call psb_error() + return + end if + return + +end subroutine psb_dbjac_aply + diff --git a/src/prec/psb_dmlprc_aply.f90 b/src/prec/psb_dmlprc_aply.f90 new file mode 100644 index 00000000..c014791d --- /dev/null +++ b/src/prec/psb_dmlprc_aply.f90 @@ -0,0 +1,529 @@ + +!!$ +!!$ +!!$ MPcube: Multilevel Parallel Preconditioners Package +!!$ for +!!$ Parallel Sparse BLAS v2.0 +!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ Daniela Di Serafino II University of Naples +!!$ Pasqua D'Ambra ICAR-CNR +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MPCUBE group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MPCUBE GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +subroutine psb_dmlprc_aply(baseprecv,x,beta,y,desc_data,trans,work,info) + ! + ! Compute Y <- beta*Y + K^-1 X + ! where K is a multilevel (actually 2-level) preconditioner stored in prec + ! + + use psb_serial_mod + use psb_descriptor_type + use psb_prec_type + use psb_psblas_mod + use psb_blacs_mod + use psb_const_mod + use psb_error_mod + implicit none + + type(psb_desc_type),intent(in) :: desc_data + type(psb_dbaseprc_type), intent(in) :: baseprecv(:) + real(kind(0.d0)),intent(in) :: beta + real(kind(0.d0)),intent(inout) :: x(:), y(:) + character :: trans + real(kind(0.d0)),target :: work(:) + integer, intent(out) :: info + + + ! Local variables + integer :: n_row,n_col + real(kind(1.d0)), allocatable :: tx(:),ty(:),t2l(:),w2l(:),& + & x2l(:),b2l(:),tz(:),tty(:) + character ::diagl, diagu + integer :: icontxt,nprow,npcol,me,mycol,i, isz, nrg,nr2l,err_act, iptype, int_err(5) + real(kind(1.d0)) :: omega + real(kind(1.d0)) :: t1, t2, t3, t4, t5, t6, t7, mpi_wtime + logical, parameter :: debug=.false., debugprt=.false. + real(kind(1.d0)), parameter :: one=1.d0, zero=0.d0 + integer :: ismth + external mpi_wtime + character(len=20) :: name, ch_err + + interface psb_baseprc_aply + subroutine psb_dbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info) + use psb_descriptor_type + use psb_prec_type + type(psb_desc_type),intent(in) :: desc_data + type(psb_dbaseprc_type), intent(in) :: prec + real(kind(0.d0)),intent(inout) :: x(:), y(:) + real(kind(0.d0)),intent(in) :: beta + character(len=1) :: trans + real(kind(0.d0)),target :: work(:) + integer, intent(out) :: info + end subroutine psb_dbaseprc_aply + end interface + + name='psb_dmlprc_aply' + info = 0 + call psb_erractionsave(err_act) + + + icontxt=desc_data%matrix_data(psb_ctxt_) + call blacs_gridinfo(icontxt,nprow,npcol,me,mycol) + + omega=baseprecv(2)%dprcparm(smooth_omega_) + ismth=baseprecv(2)%iprcparm(smth_kind_) + + select case(baseprecv(2)%iprcparm(ml_type_)) + case(no_ml_) + ! Should not really get here. + write(0,*) 'Smooth preconditioner with no multilevel in MLPRC_APLY????' + + case(add_ml_prec_) + + ! + ! Additive multilevel + ! + t1 = mpi_wtime() + n_row = desc_data%matrix_data(psb_n_row_) + n_col = baseprecv(1)%desc_data%matrix_data(psb_n_col_) + call psb_baseprc_aply(baseprecv(1),x,beta,y,desc_data,trans,work,info) + if(info /=0) goto 9999 + + nr2l = baseprecv(2)%desc_data%matrix_data(psb_n_col_) + nrg = baseprecv(2)%desc_data%matrix_data(psb_n_row_) + allocate(t2l(nr2l),w2l(nr2l),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + + t2l(:) = zero + w2l(:) = zero + + if (ismth /= no_smth_) then + ! + ! Smoothed aggregation + ! + allocate(tx(max(n_row,n_col)),ty(max(n_row,n_col)),& + & tz(max(n_row,n_col)),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + + tx(1:desc_data%matrix_data(psb_n_row_)) = x(1:desc_data%matrix_data(psb_n_row_)) + tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zero + ty(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zero + tz(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zero + + + if (baseprecv(2)%iprcparm(glb_smth_) >0) then + call psb_halo(tx,desc_data,info,work=work) + if(info /=0) goto 9999 + else + tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zero + end if + + call psb_csmm(one,baseprecv(2)%av(sm_pr_t_),tx,zero,t2l,info) + if(info /=0) goto 9999 + + else + ! + ! Raw aggregation, may take shortcut + ! + do i=1,desc_data%matrix_data(psb_n_row_) + t2l(baseprecv(2)%mlia(i)) = t2l(baseprecv(2)%mlia(i)) + x(i) + end do + + end if + + if (baseprecv(2)%iprcparm(coarse_mat_)==mat_repl_) Then + call gsum2d(icontxt,'All',t2l(1:nrg)) + else if (baseprecv(2)%iprcparm(coarse_mat_) /= mat_distr_) Then + write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',& + & baseprecv(2)%iprcparm(coarse_mat_) + endif + + w2l=t2l + call psb_baseprc_aply(baseprecv(2),w2l,zero,t2l,baseprecv(2)%desc_data,& + & 'N',work,info) + + + if (ismth /= no_smth_) then + + call psb_csmm(one,baseprecv(2)%av(sm_pr_),t2l,zero,ty,info) + if(info /=0) goto 9999 + ! + ! Finally add back into Y. + ! + call psb_axpby(one,ty,one,y,desc_data,info) + if(info /=0) goto 9999 + deallocate(tx,ty,tz) + + else + + do i=1, desc_data%matrix_data(psb_n_row_) + y(i) = y(i) + t2l(baseprecv(2)%mlia(i)) + enddo + + end if + + if (debugprt) write(0,*)' Y2: ',Y(:) + + deallocate(t2l,w2l) + + case(mult_ml_prec_) + + ! + ! Multiplicative multilevel + ! Pre/post smoothing versions. + + select case(baseprecv(2)%iprcparm(smth_pos_)) + + case(post_smooth_) + + + t1 = mpi_wtime() + n_row = desc_data%matrix_data(psb_n_row_) + n_col = baseprecv(1)%desc_data%matrix_data(psb_n_col_) + nr2l = baseprecv(2)%desc_data%matrix_data(psb_n_col_) + nrg = baseprecv(2)%desc_data%matrix_data(psb_n_row_) + allocate(t2l(nr2l),w2l(nr2l),tx(n_col),ty(n_col),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + + t2l(:) = zero + w2l(:) = zero + + ! + ! Need temp copies to handle Y<- betaY + K^-1 X + ! One of the temp copies is not strictly needed when beta==zero + ! + + if (debug) write(0,*)' mult_ml_apply omega ',omega + if (debugprt) write(0,*)' mult_ml_apply X: ',X(:) + call psb_axpby(one,x,zero,tx,desc_data,info) + if(info /=0) then + if (debug) write(0,*)' From axpby1 ',size(x),size(tx),n_row,n_col,nr2l,nrg + call psb_errpush(4010,name,a_err='axpby post_smooth 1') + goto 9999 + endif + if (ismth /= no_smth_) then + ! + ! Smoothed aggregation + ! + allocate(tz(max(n_row,n_col)),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + + + if (baseprecv(2)%iprcparm(glb_smth_) >0) then + call psb_halo(tx,desc_data,info,work=work) + if(info /=0) goto 9999 + else + tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zero + end if + + call psb_csmm(one,baseprecv(2)%av(sm_pr_t_),tx,zero,t2l,info) + if(info /=0) goto 9999 + + else + ! + ! Raw aggregation, may take shortcut + ! + do i=1,desc_data%matrix_data(psb_n_row_) + t2l(baseprecv(2)%mlia(i)) = t2l(baseprecv(2)%mlia(i)) + tx(i) + end do + end if + + if (baseprecv(2)%iprcparm(coarse_mat_)==mat_repl_) Then + call gsum2d(icontxt,'All',t2l(1:nrg)) + else if (baseprecv(2)%iprcparm(coarse_mat_) /= mat_distr_) Then + write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',& + & baseprecv(2)%iprcparm(coarse_mat_) + endif + + t6 = mpi_wtime() + w2l=t2l + call psb_baseprc_aply(baseprecv(2),w2l,zero,t2l,baseprecv(2)%desc_data,& + &'N',work,info) + if(info /=0) goto 9999 + + if (ismth /= no_smth_) then + if (ismth == smth_omg_) & + & call psb_halo(t2l,baseprecv(2)%desc_data,info,work=work) + call psb_csmm(one,baseprecv(2)%av(sm_pr_),t2l,zero,ty,info) + if(info /=0) goto 9999 + ! + ! Finally add back into Y. + ! + deallocate(tz) + else + ty(:) = zero + do i=1, desc_data%matrix_data(psb_n_row_) + ty(i) = ty(i) + t2l(baseprecv(2)%mlia(i)) + enddo + + end if + deallocate(t2l,w2l) + + call psb_spmm(-one,baseprecv(2)%aorig,ty,one,tx,desc_data,info,work=work) + if(info /=0) goto 9999 + + call psb_baseprc_aply(baseprecv(1),tx,one,ty,desc_data,trans,& + & work,info) + if(info /=0) goto 9999 + + call psb_axpby(one,ty,beta,y,desc_data,info) + if(info /=0) goto 9999 + + deallocate(tx,ty) + + + + case(pre_smooth_) + + t1 = mpi_wtime() + n_row = desc_data%matrix_data(psb_n_row_) + n_col = baseprecv(1)%desc_data%matrix_data(psb_n_col_) + nr2l = baseprecv(2)%desc_data%matrix_data(psb_n_col_) + nrg = baseprecv(2)%desc_data%matrix_data(psb_n_row_) + allocate(t2l(nr2l),w2l(nr2l),tx(n_col),ty(n_col),tty(n_col),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + + t2l(:) = zero + w2l(:) = zero + + ! + ! Need temp copies to handle Y<- betaY + K^-1 X + ! One of the temp copies is not strictly needed when beta==zero + ! + call psb_axpby(one,x,zero,tx,desc_data,info) + call psb_axpby(one,y,zero,ty,desc_data,info) + if(info /=0) goto 9999 + + call psb_baseprc_aply(baseprecv(1),x,zero,tty,desc_data,& + & trans,work,info) + if(info /=0) goto 9999 + + call psb_spmm(-one,baseprecv(2)%aorig,tty,one,tx,desc_data,info,work=work) + if(info /=0) goto 9999 + + if (ismth /= no_smth_) then + allocate(tz(max(n_row,n_col)),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + + + if (baseprecv(2)%iprcparm(glb_smth_) >0) then + call psb_halo(tx,desc_data,info,work=work) + if(info /=0) goto 9999 + else + tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zero + end if + + call psb_csmm(one,baseprecv(2)%av(sm_pr_t_),tx,zero,t2l,info) + if(info /=0) goto 9999 + + else + ! + ! Raw aggregation, may take shortcuts + ! + do i=1,desc_data%matrix_data(psb_n_row_) + t2l(baseprecv(2)%mlia(i)) = t2l(baseprecv(2)%mlia(i)) + tx(i) + end do + end if + + if (baseprecv(2)%iprcparm(coarse_mat_)==mat_repl_) Then + call gsum2d(icontxt,'All',t2l(1:nrg)) + else if (baseprecv(2)%iprcparm(coarse_mat_) /= mat_distr_) Then + write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',& + & baseprecv(2)%iprcparm(coarse_mat_) + endif + + t6 = mpi_wtime() + w2l=t2l + call psb_baseprc_aply(baseprecv(2),w2l,zero,t2l,baseprecv(2)%desc_data,& + & 'N',work,info) + if(info /=0) goto 9999 + + if (ismth /= no_smth_) then + + if (ismth == smth_omg_) & + & call psb_halo(t2l,baseprecv(2)%desc_data,info,work=work) + call psb_csmm(one,baseprecv(2)%av(sm_pr_),t2l,zero,ty,info) + if(info /=0) goto 9999 + + call psb_axpby(one,ty,one,tty,desc_data,info) + if(info /=0) goto 9999 + + deallocate(tz) + else + + do i=1, desc_data%matrix_data(psb_n_row_) + tty(i) = tty(i) + t2l(baseprecv(2)%mlia(i)) + enddo + + end if + + call psb_axpby(one,tty,beta,y,desc_data,info) + if(info /=0) goto 9999 + + deallocate(t2l,w2l,tx,ty,tty) + + + case(smooth_both_) + + t1 = mpi_wtime() + n_row = desc_data%matrix_data(psb_n_row_) + n_col = baseprecv(1)%desc_data%matrix_data(psb_n_col_) + nr2l = baseprecv(2)%desc_data%matrix_data(psb_n_col_) + nrg = baseprecv(2)%desc_data%matrix_data(psb_n_row_) + allocate(t2l(nr2l),w2l(nr2l),tx(n_col),ty(n_col),tty(n_col),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + + t2l(:) = zero + w2l(:) = zero + tx(:) = zero + ty(:) = zero + tty(:) = zero + + ! + ! Need temp copies to handle Y<- betaY + K^-1 X + ! One of the temp copies is not strictly needed when beta==zero + ! + call psb_axpby(one,x,zero,tx,desc_data,info) + call psb_axpby(one,y,zero,ty,desc_data,info) + if(info /=0) goto 9999 + + call psb_baseprc_aply(baseprecv(1),tx,zero,tty,desc_data,trans,work,info) + if(info /=0) goto 9999 + + call psb_spmm(-one,baseprecv(2)%aorig,tty,one,tx,desc_data,info,work=work) + if(info /=0) goto 9999 + + if (ismth /= no_smth_) then + if (baseprecv(2)%iprcparm(glb_smth_) >0) then + call psb_halo(tx,baseprecv(1)%desc_data,info,work=work) + if(info /=0) goto 9999 + else + tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zero + end if + + call psb_csmm(one,baseprecv(2)%av(sm_pr_t_),tx,zero,t2l,info) + if(info /=0) goto 9999 + else + ! + ! Raw aggregation, may take shortcuts + ! + do i=1,desc_data%matrix_data(psb_n_row_) + t2l(baseprecv(2)%mlia(i)) = t2l(baseprecv(2)%mlia(i)) + tx(i) + end do + end if + + + if (baseprecv(2)%iprcparm(coarse_mat_)==mat_repl_) Then + call gsum2d(icontxt,'All',t2l(1:nrg)) + else if (baseprecv(2)%iprcparm(coarse_mat_) /= mat_distr_) Then + write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',& + & baseprecv(2)%iprcparm(coarse_mat_) + endif + + + t6 = mpi_wtime() + w2l=t2l + call psb_baseprc_aply(baseprecv(2),w2l,zero,t2l,baseprecv(2)%desc_data,& + & 'N',work,info) + if(info /=0) goto 9999 + + if (ismth /= no_smth_) then + + if (ismth == smth_omg_) & + & call psb_halo(t2l,baseprecv(2)%desc_data,info,work=work) + call psb_csmm(one,baseprecv(2)%av(sm_pr_),t2l,zero,ty,info) + if(info /=0) goto 9999 + + call psb_axpby(one,ty,one,tty,desc_data,info) + if(info /=0) goto 9999 + + else + + do i=1, desc_data%matrix_data(psb_n_row_) + tty(i) = tty(i) + t2l(baseprecv(2)%mlia(i)) + enddo + + end if + + call psb_axpby(one,x,zero,tx,desc_data,info) + if(info /=0) goto 9999 + + call psb_spmm(-one,baseprecv(2)%aorig,tty,one,tx,desc_data,info,work=work) + if(info /=0) goto 9999 + call psb_baseprc_aply(baseprecv(1),tx,one,tty,desc_data,'N',work,info) + + + call psb_axpby(one,tty,beta,y,desc_data,info) + + deallocate(t2l,w2l,tx,ty,tty) + + + case default + + write(0,*) 'Unknown value for ml_smooth_pos',baseprecv(2)%iprcparm(smth_pos_) + + end select + + case default + write(0,*) me, 'Wrong mltype into PRC_APLY ',& + & baseprecv(2)%iprcparm(ml_type_) + end select + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_errpush(info,name) + call psb_erractionrestore(err_act) + if (err_act.eq.act_abort) then + call psb_error() + return + end if + return + +end subroutine psb_dmlprc_aply diff --git a/src/prec/psb_dprc_aply.f90 b/src/prec/psb_dprc_aply.f90 new file mode 100644 index 00000000..06d59c9e --- /dev/null +++ b/src/prec/psb_dprc_aply.f90 @@ -0,0 +1,261 @@ +!!$ +!!$ +!!$ MPcube: Multilevel Parallel Preconditioners Package +!!$ for +!!$ Parallel Sparse BLAS v2.0 +!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ Daniela Di Serafino II University of Naples +!!$ Pasqua D'Ambra ICAR-CNR +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MPCUBE group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MPCUBE GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +subroutine psb_dprc_aply(prec,x,y,desc_data,info,trans, work) + + use psb_serial_mod + use psb_descriptor_type + use psb_prec_type + use psb_psblas_mod + use psb_const_mod + use psb_error_mod + implicit none + + type(psb_desc_type),intent(in) :: desc_data + type(psb_dprec_type), intent(in) :: prec + real(kind(0.d0)),intent(inout) :: x(:), y(:) + integer, intent(out) :: info + character(len=1), optional :: trans + real(kind(0.d0)), optional, target :: work(:) + + ! Local variables + character ::trans_ + real(kind(1.d0)), pointer :: work_(:) + integer :: icontxt,nprow,npcol,me,mycol,err_act, int_err(5) + logical,parameter :: debug=.false., debugprt=.false. + real(kind(1.d0)), parameter :: one=1.d0, zero=0.d0 + external mpi_wtime + character(len=20) :: name, ch_err + + interface psb_baseprc_aply + subroutine psb_dbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info) + use psb_descriptor_type + use psb_prec_type + type(psb_desc_type),intent(in) :: desc_data + type(psb_dbaseprc_type), intent(in) :: prec + real(kind(0.d0)),intent(inout) :: x(:), y(:) + real(kind(0.d0)),intent(in) :: beta + character(len=1) :: trans + real(kind(0.d0)),target :: work(:) + integer, intent(out) :: info + end subroutine psb_dbaseprc_aply + end interface + + interface psb_mlprc_aply + subroutine psb_dmlprc_aply(baseprecv,x,beta,y,desc_data,trans,work,info) + use psb_descriptor_type + use psb_prec_type + type(psb_desc_type),intent(in) :: desc_data + type(psb_dbaseprc_type), intent(in) :: baseprecv(:) + real(kind(0.d0)),intent(in) :: beta + real(kind(0.d0)),intent(inout) :: x(:), y(:) + character :: trans + real(kind(0.d0)),target :: work(:) + integer, intent(out) :: info + end subroutine psb_dmlprc_aply + end interface + + name='psb_dprc_aply' + info = 0 + call psb_erractionsave(err_act) + + icontxt=desc_data%matrix_data(psb_ctxt_) + call blacs_gridinfo(icontxt,nprow,npcol,me,mycol) + + if (present(trans)) then + trans_=trans + else + trans_='N' + end if + + if (present(work)) then + work_ => work + else + allocate(work_(4*desc_data%matrix_data(psb_n_col_)),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + + end if + + if (.not.(associated(prec%baseprecv))) then + write(0,*) 'Inconsistent preconditioner: neither SMTH nor BASE?' + end if + if (size(prec%baseprecv) >1) then + if (debug) write(0,*) 'Into mlprc_aply',size(x),size(y) + call psb_mlprc_aply(prec%baseprecv,x,zero,y,desc_data,trans_,work_,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='psb_dmlprc_aply') + goto 9999 + end if + + else if (size(prec%baseprecv) == 1) then + call psb_baseprc_aply(prec%baseprecv(1),x,zero,y,desc_data,trans_, work_,info) + else + write(0,*) 'Inconsistent preconditioner: size of baseprecv???' + endif + + if (present(work)) then + else + deallocate(work_) + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act.eq.act_abort) then + call psb_error() + return + end if + return + +end subroutine psb_dprc_aply + + + +!!$ +!!$ +!!$ MPcube: Multilevel Parallel Preconditioners Package +!!$ for +!!$ Parallel Sparse BLAS v2.0 +!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ Daniela Di Serafino II University of Naples +!!$ Pasqua D'Ambra ICAR-CNR +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MPCUBE group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MPCUBE GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ + +subroutine psb_dprc_aply1(prec,x,desc_data,info,trans) + + use psb_serial_mod + use psb_descriptor_type + use psb_prec_type + use psb_psblas_mod + use psb_const_mod + use psb_error_mod + implicit none + + type(psb_desc_type),intent(in) :: desc_data + type(psb_dprec_type), intent(in) :: prec + real(kind(0.d0)),intent(inout) :: x(:) + integer, intent(out) :: info + character(len=1), optional :: trans + logical,parameter :: debug=.false., debugprt=.false. + real(kind(1.d0)), parameter :: one=1.d0, zero=0.d0 + + interface + subroutine psb_dprc_aply(prec,x,y,desc_data,info,trans, work) + + use psb_descriptor_type + use psb_prec_type + implicit none + + type(psb_desc_type),intent(in) :: desc_data + type(psb_dprec_type), intent(in) :: prec + real(kind(0.d0)),intent(inout) :: x(:), y(:) + integer, intent(out) :: info + character(len=1), optional :: trans + real(kind(0.d0)), optional, target :: work(:) + end subroutine psb_dprc_aply + end interface + + ! Local variables + character :: trans_ + integer :: icontxt,nprow,npcol,me,mycol,i, isz, err_act, int_err(5) + real(kind(1.d0)), pointer :: WW(:), w1(:) + character(len=20) :: name, ch_err + name='psb_dprec1' + info = 0 + call psb_erractionsave(err_act) + + + icontxt=desc_data%matrix_data(psb_ctxt_) + call blacs_gridinfo(icontxt,nprow,npcol,me,mycol) + if (present(trans)) then + trans_=trans + else + trans_='N' + end if + + allocate(ww(size(x)),w1(size(x)),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + if (debug) write(0,*) 'Prc_aply1 Size(x) ',size(x), size(ww),size(w1) + call psb_dprc_aply(prec,x,ww,desc_data,info,trans_,work=w1) + if(info /=0) goto 9999 + x(:) = ww(:) + deallocate(ww,W1) + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_errpush(info,name) + call psb_erractionrestore(err_act) + if (err_act.eq.act_abort) then + call psb_error() + return + end if + return +end subroutine psb_dprc_aply1 diff --git a/src/prec/psb_dprec.f90 b/src/prec/psb_dprec.f90 deleted file mode 100644 index 66eaf277..00000000 --- a/src/prec/psb_dprec.f90 +++ /dev/null @@ -1,1384 +0,0 @@ -!!$ -!!$ -!!$ MPcube: Multilevel Parallel Preconditioners Package -!!$ for -!!$ Parallel Sparse BLAS v2.0 -!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata -!!$ Alfredo Buttari University of Rome Tor Vergata -!!$ Daniela Di Serafino II University of Naples -!!$ Pasqua D'Ambra ICAR-CNR -!!$ -!!$ Redistribution and use in source and binary forms, with or without -!!$ modification, are permitted provided that the following conditions -!!$ are met: -!!$ 1. Redistributions of source code must retain the above copyright -!!$ notice, this list of conditions and the following disclaimer. -!!$ 2. Redistributions in binary form must reproduce the above copyright -!!$ notice, this list of conditions, and the following disclaimer in the -!!$ documentation and/or other materials provided with the distribution. -!!$ 3. The name of the MPCUBE group or the names of its contributors may -!!$ not be used to endorse or promote products derived from this -!!$ software without specific written permission. -!!$ -!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR -!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MPCUBE GROUP OR ITS CONTRIBUTORS -!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR -!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF -!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS -!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN -!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -!!$ POSSIBILITY OF SUCH DAMAGE. -!!$ -!!$ -subroutine psb_dprc_aply(prec,x,y,desc_data,info,trans, work) - - use psb_serial_mod - use psb_descriptor_type - use psb_prec_type - use psb_psblas_mod - use psb_const_mod - use psb_error_mod - implicit none - - type(psb_desc_type),intent(in) :: desc_data - type(psb_dprec_type), intent(in) :: prec - real(kind(0.d0)),intent(inout) :: x(:), y(:) - integer, intent(out) :: info - character(len=1), optional :: trans - real(kind(0.d0)), optional, target :: work(:) - - ! Local variables - character ::trans_ - real(kind(1.d0)), pointer :: work_(:) - integer :: icontxt,nprow,npcol,me,mycol,err_act, int_err(5) - logical,parameter :: debug=.false., debugprt=.false. - real(kind(1.d0)), parameter :: one=1.d0, zero=0.d0 - external mpi_wtime - character(len=20) :: name, ch_err - - interface psb_baseprc_aply - subroutine psb_dbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info) - use psb_descriptor_type - use psb_prec_type - type(psb_desc_type),intent(in) :: desc_data - type(psb_dbaseprc_type), intent(in) :: prec - real(kind(0.d0)),intent(inout) :: x(:), y(:) - real(kind(0.d0)),intent(in) :: beta - character(len=1) :: trans - real(kind(0.d0)),target :: work(:) - integer, intent(out) :: info - end subroutine psb_dbaseprc_aply - end interface - - interface psb_mlprc_aply - subroutine psb_dmlprc_aply(baseprecv,x,beta,y,desc_data,trans,work,info) - use psb_descriptor_type - use psb_prec_type - type(psb_desc_type),intent(in) :: desc_data - type(psb_dbaseprc_type), intent(in) :: baseprecv(:) - real(kind(0.d0)),intent(in) :: beta - real(kind(0.d0)),intent(inout) :: x(:), y(:) - character :: trans - real(kind(0.d0)),target :: work(:) - integer, intent(out) :: info - end subroutine psb_dmlprc_aply - end interface - - name='psb_dprc_aply' - info = 0 - call psb_erractionsave(err_act) - - icontxt=desc_data%matrix_data(psb_ctxt_) - call blacs_gridinfo(icontxt,nprow,npcol,me,mycol) - - if (present(trans)) then - trans_=trans - else - trans_='N' - end if - - if (present(work)) then - work_ => work - else - allocate(work_(4*desc_data%matrix_data(psb_n_col_)),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if - - end if - - if (.not.(associated(prec%baseprecv))) then - write(0,*) 'Inconsistent preconditioner: neither SMTH nor BASE?' - end if - if (size(prec%baseprecv) >1) then - if (debug) write(0,*) 'Into mlprc_aply',size(x),size(y) - call psb_mlprc_aply(prec%baseprecv,x,zero,y,desc_data,trans_,work_,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_dmlprc_aply') - goto 9999 - end if - - else if (size(prec%baseprecv) == 1) then - call psb_baseprc_aply(prec%baseprecv(1),x,zero,y,desc_data,trans_, work_,info) - else - write(0,*) 'Inconsistent preconditioner: size of baseprecv???' - endif - - if (present(work)) then - else - deallocate(work_) - end if - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act.eq.act_abort) then - call psb_error() - return - end if - return - -end subroutine psb_dprc_aply - - -!!$ -!!$ -!!$ MPcube: Multilevel Parallel Preconditioners Package -!!$ for -!!$ Parallel Sparse BLAS v2.0 -!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata -!!$ Alfredo Buttari University of Rome Tor Vergata -!!$ Daniela Di Serafino II University of Naples -!!$ Pasqua D'Ambra ICAR-CNR -!!$ -!!$ Redistribution and use in source and binary forms, with or without -!!$ modification, are permitted provided that the following conditions -!!$ are met: -!!$ 1. Redistributions of source code must retain the above copyright -!!$ notice, this list of conditions and the following disclaimer. -!!$ 2. Redistributions in binary form must reproduce the above copyright -!!$ notice, this list of conditions, and the following disclaimer in the -!!$ documentation and/or other materials provided with the distribution. -!!$ 3. The name of the MPCUBE group or the names of its contributors may -!!$ not be used to endorse or promote products derived from this -!!$ software without specific written permission. -!!$ -!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR -!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MPCUBE GROUP OR ITS CONTRIBUTORS -!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR -!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF -!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS -!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN -!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -!!$ POSSIBILITY OF SUCH DAMAGE. -!!$ -!!$ - -subroutine psb_dbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info) - ! - ! Compute Y <- beta*Y + K^-1 X - ! where K is a a basic preconditioner stored in prec - ! - - use psb_serial_mod - use psb_descriptor_type - use psb_prec_type - use psb_psblas_mod - use psb_const_mod - use psb_error_mod - implicit none - - type(psb_desc_type),intent(in) :: desc_data - type(psb_dbaseprc_type), intent(in) :: prec - real(kind(0.d0)),intent(inout) :: x(:), y(:) - real(kind(0.d0)),intent(in) :: beta - character(len=1) :: trans - real(kind(0.d0)),target :: work(:) - integer, intent(out) :: info - - ! Local variables - integer :: n_row,n_col, int_err(5) - real(kind(1.d0)), pointer :: ww(:), aux(:), tx(:),ty(:) - character ::diagl, diagu - integer :: icontxt,nprow,npcol,me,mycol,i, isz, nrg, err_act - real(kind(1.d0)) :: t1, t2, t3, t4, t5, t6, t7, mpi_wtime - logical,parameter :: debug=.false., debugprt=.false. - real(kind(1.d0)), parameter :: one=1.d0, zero=0.d0 - external mpi_wtime - character(len=20) :: name, ch_err - - interface psb_bjac_aply - subroutine psb_dbjac_aply(prec,x,beta,y,desc_data,trans,work,info) - use psb_descriptor_type - use psb_prec_type - type(psb_desc_type), intent(in) :: desc_data - type(psb_dbaseprc_type), intent(in) :: prec - real(kind(0.d0)),intent(inout) :: x(:), y(:) - real(kind(0.d0)),intent(in) :: beta - character(len=1) :: trans - real(kind(0.d0)),target :: work(:) - integer, intent(out) :: info - end subroutine psb_dbjac_aply - end interface - - name='psb_dbaseprc_aply' - info = 0 - call psb_erractionsave(err_act) - - icontxt=desc_data%matrix_data(psb_ctxt_) - call blacs_gridinfo(icontxt,nprow,npcol,me,mycol) - - diagl='U' - diagu='U' - - select case(trans) - case('N','n') - case('T','t','C','c') - case default - info=40 - int_err(1)=6 - ch_err(2:2)=trans - goto 9999 - end select - - select case(prec%iprcparm(p_type_)) - - case(noprec_) - - n_row=desc_data%matrix_data(psb_n_row_) - if (beta==zero) then - y(1:n_row) = x(1:n_row) - else if (beta==one) then - y(1:n_row) = x(1:n_row) + y(1:n_row) - else if (beta==-one) then - y(1:n_row) = x(1:n_row) - y(1:n_row) - else - y(1:n_row) = x(1:n_row) + beta*y(1:n_row) - end if - - case(diagsc_) - - n_row=desc_data%matrix_data(psb_n_row_) - if (beta==zero) then - y(1:n_row) = x(1:n_row)*prec%d(1:n_row) - else if (beta==one) then - y(1:n_row) = x(1:n_row)*prec%d(1:n_row) + y(1:n_row) - else if (beta==-one) then - y(1:n_row) = x(1:n_row)*prec%d(1:n_row) - y(1:n_row) - else - y(1:n_row) = x(1:n_row)*prec%d(1:n_row) + beta*y(1:n_row) - end if - - case(bja_) - - call psb_bjac_aply(prec,x,beta,y,desc_data,trans,work,info) - if(info.ne.0) then - info=4010 - ch_err='psb_bjac_aply' - goto 9999 - end if - - case(asm_,ras_,ash_,rash_) - - if (prec%iprcparm(n_ovr_)==0) then - ! shortcut: this fixes performance for RAS(0) == BJA - call psb_bjac_aply(prec,x,beta,y,desc_data,trans,work,info) - if(info.ne.0) then - info=4010 - ch_err='psb_bjacaply' - goto 9999 - end if - - else - ! Note: currently trans is unused. - n_row=prec%desc_data%matrix_data(psb_n_row_) - n_col=prec%desc_data%matrix_data(psb_n_col_) - - isz=max(n_row,N_COL) - if ((6*isz) <= size(work)) then - ww => work(1:isz) - tx => work(isz+1:2*isz) - ty => work(2*isz+1:3*isz) - aux => work(3*isz+1:) - else if ((4*isz) <= size(work)) then - aux => work(1:) - allocate(ww(isz),tx(isz),ty(isz),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if - else if ((3*isz) <= size(work)) then - ww => work(1:isz) - tx => work(isz+1:2*isz) - ty => work(2*isz+1:3*isz) - allocate(aux(4*isz),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if - - else - allocate(ww(isz),tx(isz),ty(isz),& - &aux(4*isz),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if - - endif - - if (debugprt) write(0,*)' vdiag: ',prec%d(:) - if (debug) write(0,*) 'Bi-CGSTAB with Additive Schwarz prec' - - tx(1:desc_data%matrix_data(psb_n_row_)) = x(1:desc_data%matrix_data(psb_n_row_)) - tx(desc_data%matrix_data(psb_n_row_)+1:isz) = zero - - if (prec%iprcparm(restr_)==psb_halo_) then - call psb_halo(tx,prec%desc_data,info,work=aux) - if(info /=0) then - info=4010 - ch_err='psb_halo' - goto 9999 - end if - else if (prec%iprcparm(restr_) /= psb_none_) then - write(0,*) 'Problem in PRC_APLY: Unknown value for restriction ',& - &prec%iprcparm(restr_) - end if - - if (prec%iprcparm(iren_)>0) then - call dgelp('N',n_row,1,prec%perm,tx,isz,ww,isz,info) - if(info /=0) then - info=4010 - ch_err='psb_dgelp' - goto 9999 - end if - endif - - call psb_bjac_aply(prec,tx,zero,ty,prec%desc_data,trans,aux,info) - if(info.ne.0) then - info=4010 - ch_err='psb_bjac_aply' - goto 9999 - end if - - if (prec%iprcparm(iren_)>0) then - call dgelp('N',n_row,1,prec%invperm,ty,isz,ww,isz,info) - if(info /=0) then - info=4010 - ch_err='psb_dgelp' - goto 9999 - end if - endif - - select case (prec%iprcparm(prol_)) - - case(psb_none_) - ! Would work anyway, but since it's supposed to do nothing... - ! call f90_psovrl(ty,prec%desc_data,update_type=prec%a_restrict) - - case(psb_sum_,psb_avg_) - call psb_ovrl(ty,prec%desc_data,info,& - & update_type=prec%iprcparm(prol_),work=aux) - if(info /=0) then - info=4010 - ch_err='psb_ovrl' - goto 9999 - end if - - case default - write(0,*) 'Problem in PRC_APLY: Unknown value for prolongation ',& - & prec%iprcparm(prol_) - end select - - if (beta == zero) then - y(1:desc_data%matrix_data(psb_n_row_)) = ty(1:desc_data%matrix_data(psb_n_row_)) - else if (beta == one) then - y(1:desc_data%matrix_data(psb_n_row_)) = y(1:desc_data%matrix_data(psb_n_row_)) + & - & ty(1:desc_data%matrix_data(psb_n_row_)) - else if (beta == -one) then - y(1:desc_data%matrix_data(psb_n_row_)) = -y(1:desc_data%matrix_data(psb_n_row_)) + & - & ty(1:desc_data%matrix_data(psb_n_row_)) - else - y(1:desc_data%matrix_data(psb_n_row_)) = beta*y(1:desc_data%matrix_data(psb_n_row_)) + & - & ty(1:desc_data%matrix_data(psb_n_row_)) - end if - - - if ((6*isz) <= size(work)) then - else if ((4*isz) <= size(work)) then - deallocate(ww,tx,ty) - else if ((3*isz) <= size(work)) then - deallocate(aux) - else - deallocate(ww,aux,tx,ty) - endif - end if - case default - write(0,*) 'Invalid PRE%PREC ',prec%iprcparm(p_type_),':',& - & min_prec_,noprec_,diagsc_,bja_,& - & asm_,ras_,ash_,rash_ - end select - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_errpush(info,name,i_err=int_err,a_err=ch_err) - call psb_erractionrestore(err_act) - if (err_act.eq.act_abort) then - call psb_error() - return - end if - return - -end subroutine psb_dbaseprc_aply - - - -!!$ -!!$ -!!$ MPcube: Multilevel Parallel Preconditioners Package -!!$ for -!!$ Parallel Sparse BLAS v2.0 -!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata -!!$ Alfredo Buttari University of Rome Tor Vergata -!!$ Daniela Di Serafino II University of Naples -!!$ Pasqua D'Ambra ICAR-CNR -!!$ -!!$ Redistribution and use in source and binary forms, with or without -!!$ modification, are permitted provided that the following conditions -!!$ are met: -!!$ 1. Redistributions of source code must retain the above copyright -!!$ notice, this list of conditions and the following disclaimer. -!!$ 2. Redistributions in binary form must reproduce the above copyright -!!$ notice, this list of conditions, and the following disclaimer in the -!!$ documentation and/or other materials provided with the distribution. -!!$ 3. The name of the MPCUBE group or the names of its contributors may -!!$ not be used to endorse or promote products derived from this -!!$ software without specific written permission. -!!$ -!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR -!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MPCUBE GROUP OR ITS CONTRIBUTORS -!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR -!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF -!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS -!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN -!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -!!$ POSSIBILITY OF SUCH DAMAGE. -!!$ -!!$ -subroutine psb_dbjac_aply(prec,x,beta,y,desc_data,trans,work,info) - ! - ! Compute Y <- beta*Y + K^-1 X - ! where K is a a Block Jacobi preconditioner stored in prec - ! Note that desc_data may or may not be the same as prec%desc_data, - ! but since both are INTENT(IN) this should be legal. - ! - - use psb_serial_mod - use psb_descriptor_type - use psb_prec_type - use psb_psblas_mod - use psb_const_mod - use psb_error_mod - implicit none - - type(psb_desc_type), intent(in) :: desc_data - type(psb_dbaseprc_type), intent(in) :: prec - real(kind(0.d0)),intent(inout) :: x(:), y(:) - real(kind(0.d0)),intent(in) :: beta - character(len=1) :: trans - real(kind(0.d0)),target :: work(:) - integer, intent(out) :: info - - ! Local variables - integer :: n_row,n_col - real(kind(1.d0)), pointer :: ww(:), aux(:), tx(:),ty(:),tb(:) - character ::diagl, diagu - integer :: icontxt,nprow,npcol,me,mycol,i, isz, nrg, err_act, int_err(5) - real(kind(1.d0)) :: t1, t2, t3, t4, t5, t6, t7, mpi_wtime - logical,parameter :: debug=.false., debugprt=.false. - real(kind(1.d0)), parameter :: one=1.d0, zero=0.d0 - external mpi_wtime - character(len=20) :: name, ch_err - - name='psb_bjac_aply' - info = 0 - call psb_erractionsave(err_act) - - icontxt=desc_data%matrix_data(psb_ctxt_) - call blacs_gridinfo(icontxt,nprow,npcol,me,mycol) - - diagl='U' - diagu='U' - - select case(trans) - case('N','n') - case('T','t','C','c') - case default - call psb_errpush(40,name) - goto 9999 - end select - - - n_row=desc_data%matrix_data(psb_n_row_) - n_col=desc_data%matrix_data(psb_n_col_) - - if (n_col <= size(work)) then - ww => work(1:n_col) - if ((4*n_col+n_col) <= size(work)) then - aux => work(n_col+1:) - else - allocate(aux(4*n_col),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if - - endif - else - allocate(ww(n_col),aux(4*n_col),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if - endif - - - if (prec%iprcparm(jac_sweeps_) == 1) then - - - select case(prec%iprcparm(f_type_)) - case(f_ilu_n_,f_ilu_e_) - - select case(trans) - case('N','n') - - call psb_spsm(one,prec%av(l_pr_),x,zero,ww,desc_data,info,& - & trans='N',unit=diagl,choice=psb_none_,work=aux) - if(info /=0) goto 9999 - ww(1:n_row) = ww(1:n_row)*prec%d(1:n_row) - call psb_spsm(one,prec%av(u_pr_),ww,beta,y,desc_data,info,& - & trans='N',unit=diagu,choice=psb_none_, work=aux) - if(info /=0) goto 9999 - - case('T','t','C','c') - call psb_spsm(one,prec%av(u_pr_),x,zero,ww,desc_data,info,& - & trans=trans,unit=diagu,choice=psb_none_, work=aux) - if(info /=0) goto 9999 - ww(1:n_row) = ww(1:n_row)*prec%d(1:n_row) - call psb_spsm(one,prec%av(l_pr_),ww,beta,y,desc_data,info,& - & trans=trans,unit=diagl,choice=psb_none_,work=aux) - if(info /=0) goto 9999 - - end select - - case(f_slu_) - - ww(1:n_row) = x(1:n_row) - - select case(trans) - case('N','n') - call psb_slu_solve(0,n_row,1,ww,n_row,prec%iprcparm(slu_ptr_),info) - case('T','t','C','c') - call psb_slu_solve(1,n_row,1,ww,n_row,prec%iprcparm(slu_ptr_),info) - end select - - if(info /=0) goto 9999 - - if (beta == 0.d0) then - y(1:n_row) = ww(1:n_row) - else if (beta==1.d0) then - y(1:n_row) = ww(1:n_row) + y(1:n_row) - else if (beta==-1.d0) then - y(1:n_row) = ww(1:n_row) - y(1:n_row) - else - y(1:n_row) = ww(1:n_row) + beta*y(1:n_row) - endif - case (f_umf_) - - - select case(trans) - case('N','n') - call psb_umf_solve(0,n_row,ww,x,n_row,prec%iprcparm(umf_numptr_),info) - case('T','t','C','c') - call psb_umf_solve(1,n_row,ww,x,n_row,prec%iprcparm(umf_numptr_),info) - end select - - if(info /=0) goto 9999 - - if (beta == 0.d0) then - y(1:n_row) = ww(1:n_row) - else if (beta==1.d0) then - y(1:n_row) = ww(1:n_row) + y(1:n_row) - else if (beta==-1.d0) then - y(1:n_row) = ww(1:n_row) - y(1:n_row) - else - y(1:n_row) = ww(1:n_row) + beta*y(1:n_row) - endif - - case default - write(0,*) 'Unknown factorization type in bjac_aply',prec%iprcparm(f_type_) - end select - if (debugprt) write(0,*)' Y: ',y(:) - - else if (prec%iprcparm(jac_sweeps_) > 1) then - - ! Note: we have to add TRANS to this one !!!!!!!!! - - if (size(prec%av) < ap_nd_) then - info = 4011 - goto 9999 - endif - - allocate(tx(n_col),ty(n_col),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if - - tx = zero - ty = zero - select case(prec%iprcparm(f_type_)) - case(f_ilu_n_,f_ilu_e_) - do i=1, prec%iprcparm(jac_sweeps_) - ! X(k+1) = M^-1*(b-N*X(k)) - ty(1:n_row) = x(1:n_row) - call psb_spmm(-one,prec%av(ap_nd_),tx,one,ty,& - & prec%desc_data,info,work=aux) - if(info /=0) goto 9999 - call psb_spsm(one,prec%av(l_pr_),ty,zero,ww,& - & prec%desc_data,info,& - & trans='N',unit='U',choice=psb_none_,work=aux) - if(info /=0) goto 9999 - ww(1:n_row) = ww(1:n_row)*prec%d(1:n_row) - call psb_spsm(one,prec%av(u_pr_),ww,zero,tx,& - & prec%desc_data,info,& - & trans='N',unit='U',choice=psb_none_,work=aux) - if(info /=0) goto 9999 - end do - - case(f_slu_) - do i=1, prec%iprcparm(jac_sweeps_) - ! X(k+1) = M^-1*(b-N*X(k)) - ty(1:n_row) = x(1:n_row) - call psb_spmm(-one,prec%av(ap_nd_),tx,one,ty,& - & prec%desc_data,info,work=aux) - if(info /=0) goto 9999 - - call psb_slu_solve(0,n_row,1,ty,n_row,prec%iprcparm(slu_ptr_),info) - if(info /=0) goto 9999 - tx(1:n_row) = ty(1:n_row) - end do - case(f_umf_) - do i=1, prec%iprcparm(jac_sweeps_) - ! X(k+1) = M^-1*(b-N*X(k)) - ty(1:n_row) = x(1:n_row) - call psb_spmm(-one,prec%av(ap_nd_),tx,one,ty,& - & prec%desc_data,info,work=aux) - if(info /=0) goto 9999 - - call psb_umf_solve(0,n_row,ww,ty,n_row,& - & prec%iprcparm(umf_numptr_),info) - if(info /=0) goto 9999 - tx(1:n_row) = ww(1:n_row) - end do - - end select - - if (beta == 0.d0) then - y(1:n_row) = tx(1:n_row) - else if (beta==1.d0) then - y(1:n_row) = tx(1:n_row) + y(1:n_row) - else if (beta==-1.d0) then - y(1:n_row) = tx(1:n_row) - y(1:n_row) - else - y(1:n_row) = tx(1:n_row) + beta*y(1:n_row) - endif - - deallocate(tx,ty) - - - else - - goto 9999 - - endif - - if (n_col <= size(work)) then - if ((4*n_col+n_col) <= size(work)) then - else - deallocate(aux) - endif - else - deallocate(ww,aux) - endif - - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_errpush(info,name,i_err=int_err,a_err=ch_err) - call psb_erractionrestore(err_act) - if (err_act.eq.act_abort) then - call psb_error() - return - end if - return - -end subroutine psb_dbjac_aply - - -!!$ -!!$ -!!$ MPcube: Multilevel Parallel Preconditioners Package -!!$ for -!!$ Parallel Sparse BLAS v2.0 -!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata -!!$ Alfredo Buttari University of Rome Tor Vergata -!!$ Daniela Di Serafino II University of Naples -!!$ Pasqua D'Ambra ICAR-CNR -!!$ -!!$ Redistribution and use in source and binary forms, with or without -!!$ modification, are permitted provided that the following conditions -!!$ are met: -!!$ 1. Redistributions of source code must retain the above copyright -!!$ notice, this list of conditions and the following disclaimer. -!!$ 2. Redistributions in binary form must reproduce the above copyright -!!$ notice, this list of conditions, and the following disclaimer in the -!!$ documentation and/or other materials provided with the distribution. -!!$ 3. The name of the MPCUBE group or the names of its contributors may -!!$ not be used to endorse or promote products derived from this -!!$ software without specific written permission. -!!$ -!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR -!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MPCUBE GROUP OR ITS CONTRIBUTORS -!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR -!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF -!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS -!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN -!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -!!$ POSSIBILITY OF SUCH DAMAGE. -!!$ -!!$ -subroutine psb_dmlprc_aply(baseprecv,x,beta,y,desc_data,trans,work,info) - ! - ! Compute Y <- beta*Y + K^-1 X - ! where K is a multilevel (actually 2-level) preconditioner stored in prec - ! - - use psb_serial_mod - use psb_descriptor_type - use psb_prec_type - use psb_psblas_mod - use psb_blacs_mod - use psb_const_mod - use psb_error_mod - implicit none - - type(psb_desc_type),intent(in) :: desc_data - type(psb_dbaseprc_type), intent(in) :: baseprecv(:) - real(kind(0.d0)),intent(in) :: beta - real(kind(0.d0)),intent(inout) :: x(:), y(:) - character :: trans - real(kind(0.d0)),target :: work(:) - integer, intent(out) :: info - - - ! Local variables - integer :: n_row,n_col - real(kind(1.d0)), allocatable :: tx(:),ty(:),t2l(:),w2l(:),& - & x2l(:),b2l(:),tz(:),tty(:) - character ::diagl, diagu - integer :: icontxt,nprow,npcol,me,mycol,i, isz, nrg,nr2l,err_act, iptype, int_err(5) - real(kind(1.d0)) :: omega - real(kind(1.d0)) :: t1, t2, t3, t4, t5, t6, t7, mpi_wtime - logical, parameter :: debug=.false., debugprt=.false. - real(kind(1.d0)), parameter :: one=1.d0, zero=0.d0 - integer :: ismth - external mpi_wtime - character(len=20) :: name, ch_err - - interface psb_baseprc_aply - subroutine psb_dbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info) - use psb_descriptor_type - use psb_prec_type - type(psb_desc_type),intent(in) :: desc_data - type(psb_dbaseprc_type), intent(in) :: prec - real(kind(0.d0)),intent(inout) :: x(:), y(:) - real(kind(0.d0)),intent(in) :: beta - character(len=1) :: trans - real(kind(0.d0)),target :: work(:) - integer, intent(out) :: info - end subroutine psb_dbaseprc_aply - end interface - - name='psb_dmlprc_aply' - info = 0 - call psb_erractionsave(err_act) - - - icontxt=desc_data%matrix_data(psb_ctxt_) - call blacs_gridinfo(icontxt,nprow,npcol,me,mycol) - - omega=baseprecv(2)%dprcparm(smooth_omega_) - ismth=baseprecv(2)%iprcparm(smth_kind_) - - select case(baseprecv(2)%iprcparm(ml_type_)) - case(no_ml_) - ! Should not really get here. - write(0,*) 'Smooth preconditioner with no multilevel in MLPRC_APLY????' - - case(add_ml_prec_) - - ! - ! Additive multilevel - ! - t1 = mpi_wtime() - n_row = desc_data%matrix_data(psb_n_row_) - n_col = baseprecv(1)%desc_data%matrix_data(psb_n_col_) - call psb_baseprc_aply(baseprecv(1),x,beta,y,desc_data,trans,work,info) - if(info /=0) goto 9999 - - nr2l = baseprecv(2)%desc_data%matrix_data(psb_n_col_) - nrg = baseprecv(2)%desc_data%matrix_data(psb_n_row_) - allocate(t2l(nr2l),w2l(nr2l),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if - - t2l(:) = zero - w2l(:) = zero - - if (ismth /= no_smth_) then - ! - ! Smoothed aggregation - ! - allocate(tx(max(n_row,n_col)),ty(max(n_row,n_col)),& - & tz(max(n_row,n_col)),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if - - tx(1:desc_data%matrix_data(psb_n_row_)) = x(1:desc_data%matrix_data(psb_n_row_)) - tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zero - ty(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zero - tz(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zero - - - if (baseprecv(2)%iprcparm(glb_smth_) >0) then - call psb_halo(tx,desc_data,info,work=work) - if(info /=0) goto 9999 - else - tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zero - end if - - call psb_csmm(one,baseprecv(2)%av(sm_pr_t_),tx,zero,t2l,info) - if(info /=0) goto 9999 - - else - ! - ! Raw aggregation, may take shortcut - ! - do i=1,desc_data%matrix_data(psb_n_row_) - t2l(baseprecv(2)%mlia(i)) = t2l(baseprecv(2)%mlia(i)) + x(i) - end do - - end if - - if (baseprecv(2)%iprcparm(coarse_mat_)==mat_repl_) Then - call gsum2d(icontxt,'All',t2l(1:nrg)) - else if (baseprecv(2)%iprcparm(coarse_mat_) /= mat_distr_) Then - write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',& - & baseprecv(2)%iprcparm(coarse_mat_) - endif - - w2l=t2l - call psb_baseprc_aply(baseprecv(2),w2l,zero,t2l,baseprecv(2)%desc_data,& - & 'N',work,info) - - - if (ismth /= no_smth_) then - - call psb_csmm(one,baseprecv(2)%av(sm_pr_),t2l,zero,ty,info) - if(info /=0) goto 9999 - ! - ! Finally add back into Y. - ! - call psb_axpby(one,ty,one,y,desc_data,info) - if(info /=0) goto 9999 - deallocate(tx,ty,tz) - - else - - do i=1, desc_data%matrix_data(psb_n_row_) - y(i) = y(i) + t2l(baseprecv(2)%mlia(i)) - enddo - - end if - - if (debugprt) write(0,*)' Y2: ',Y(:) - - deallocate(t2l,w2l) - - case(mult_ml_prec_) - - ! - ! Multiplicative multilevel - ! Pre/post smoothing versions. - - select case(baseprecv(2)%iprcparm(smth_pos_)) - - case(post_smooth_) - - - t1 = mpi_wtime() - n_row = desc_data%matrix_data(psb_n_row_) - n_col = baseprecv(1)%desc_data%matrix_data(psb_n_col_) - nr2l = baseprecv(2)%desc_data%matrix_data(psb_n_col_) - nrg = baseprecv(2)%desc_data%matrix_data(psb_n_row_) - allocate(t2l(nr2l),w2l(nr2l),tx(n_col),ty(n_col),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if - - t2l(:) = zero - w2l(:) = zero - - ! - ! Need temp copies to handle Y<- betaY + K^-1 X - ! One of the temp copies is not strictly needed when beta==zero - ! - - if (debug) write(0,*)' mult_ml_apply omega ',omega - if (debugprt) write(0,*)' mult_ml_apply X: ',X(:) - call psb_axpby(one,x,zero,tx,desc_data,info) - if(info /=0) then - if (debug) write(0,*)' From axpby1 ',size(x),size(tx),n_row,n_col,nr2l,nrg - call psb_errpush(4010,name,a_err='axpby post_smooth 1') - goto 9999 - endif - if (ismth /= no_smth_) then - ! - ! Smoothed aggregation - ! - allocate(tz(max(n_row,n_col)),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if - - - if (baseprecv(2)%iprcparm(glb_smth_) >0) then - call psb_halo(tx,desc_data,info,work=work) - if(info /=0) goto 9999 - else - tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zero - end if - - call psb_csmm(one,baseprecv(2)%av(sm_pr_t_),tx,zero,t2l,info) - if(info /=0) goto 9999 - - else - ! - ! Raw aggregation, may take shortcut - ! - do i=1,desc_data%matrix_data(psb_n_row_) - t2l(baseprecv(2)%mlia(i)) = t2l(baseprecv(2)%mlia(i)) + tx(i) - end do - end if - - if (baseprecv(2)%iprcparm(coarse_mat_)==mat_repl_) Then - call gsum2d(icontxt,'All',t2l(1:nrg)) - else if (baseprecv(2)%iprcparm(coarse_mat_) /= mat_distr_) Then - write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',& - & baseprecv(2)%iprcparm(coarse_mat_) - endif - - t6 = mpi_wtime() - w2l=t2l - call psb_baseprc_aply(baseprecv(2),w2l,zero,t2l,baseprecv(2)%desc_data,& - &'N',work,info) - if(info /=0) goto 9999 - - if (ismth /= no_smth_) then - if (ismth == smth_omg_) & - & call psb_halo(t2l,baseprecv(2)%desc_data,info,work=work) - call psb_csmm(one,baseprecv(2)%av(sm_pr_),t2l,zero,ty,info) - if(info /=0) goto 9999 - ! - ! Finally add back into Y. - ! - deallocate(tz) - else - ty(:) = zero - do i=1, desc_data%matrix_data(psb_n_row_) - ty(i) = ty(i) + t2l(baseprecv(2)%mlia(i)) - enddo - - end if - deallocate(t2l,w2l) - - call psb_spmm(-one,baseprecv(2)%aorig,ty,one,tx,desc_data,info,work=work) - if(info /=0) goto 9999 - - call psb_baseprc_aply(baseprecv(1),tx,one,ty,desc_data,trans,& - & work,info) - if(info /=0) goto 9999 - - call psb_axpby(one,ty,beta,y,desc_data,info) - if(info /=0) goto 9999 - - deallocate(tx,ty) - - - - case(pre_smooth_) - - t1 = mpi_wtime() - n_row = desc_data%matrix_data(psb_n_row_) - n_col = baseprecv(1)%desc_data%matrix_data(psb_n_col_) - nr2l = baseprecv(2)%desc_data%matrix_data(psb_n_col_) - nrg = baseprecv(2)%desc_data%matrix_data(psb_n_row_) - allocate(t2l(nr2l),w2l(nr2l),tx(n_col),ty(n_col),tty(n_col),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if - - t2l(:) = zero - w2l(:) = zero - - ! - ! Need temp copies to handle Y<- betaY + K^-1 X - ! One of the temp copies is not strictly needed when beta==zero - ! - call psb_axpby(one,x,zero,tx,desc_data,info) - call psb_axpby(one,y,zero,ty,desc_data,info) - if(info /=0) goto 9999 - - call psb_baseprc_aply(baseprecv(1),x,zero,tty,desc_data,& - & trans,work,info) - if(info /=0) goto 9999 - - call psb_spmm(-one,baseprecv(2)%aorig,tty,one,tx,desc_data,info,work=work) - if(info /=0) goto 9999 - - if (ismth /= no_smth_) then - allocate(tz(max(n_row,n_col)),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if - - - if (baseprecv(2)%iprcparm(glb_smth_) >0) then - call psb_halo(tx,desc_data,info,work=work) - if(info /=0) goto 9999 - else - tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zero - end if - - call psb_csmm(one,baseprecv(2)%av(sm_pr_t_),tx,zero,t2l,info) - if(info /=0) goto 9999 - - else - ! - ! Raw aggregation, may take shortcuts - ! - do i=1,desc_data%matrix_data(psb_n_row_) - t2l(baseprecv(2)%mlia(i)) = t2l(baseprecv(2)%mlia(i)) + tx(i) - end do - end if - - if (baseprecv(2)%iprcparm(coarse_mat_)==mat_repl_) Then - call gsum2d(icontxt,'All',t2l(1:nrg)) - else if (baseprecv(2)%iprcparm(coarse_mat_) /= mat_distr_) Then - write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',& - & baseprecv(2)%iprcparm(coarse_mat_) - endif - - t6 = mpi_wtime() - w2l=t2l - call psb_baseprc_aply(baseprecv(2),w2l,zero,t2l,baseprecv(2)%desc_data,& - & 'N',work,info) - if(info /=0) goto 9999 - - if (ismth /= no_smth_) then - - if (ismth == smth_omg_) & - & call psb_halo(t2l,baseprecv(2)%desc_data,info,work=work) - call psb_csmm(one,baseprecv(2)%av(sm_pr_),t2l,zero,ty,info) - if(info /=0) goto 9999 - - call psb_axpby(one,ty,one,tty,desc_data,info) - if(info /=0) goto 9999 - - deallocate(tz) - else - - do i=1, desc_data%matrix_data(psb_n_row_) - tty(i) = tty(i) + t2l(baseprecv(2)%mlia(i)) - enddo - - end if - - call psb_axpby(one,tty,beta,y,desc_data,info) - if(info /=0) goto 9999 - - deallocate(t2l,w2l,tx,ty,tty) - - - case(smooth_both_) - - t1 = mpi_wtime() - n_row = desc_data%matrix_data(psb_n_row_) - n_col = baseprecv(1)%desc_data%matrix_data(psb_n_col_) - nr2l = baseprecv(2)%desc_data%matrix_data(psb_n_col_) - nrg = baseprecv(2)%desc_data%matrix_data(psb_n_row_) - allocate(t2l(nr2l),w2l(nr2l),tx(n_col),ty(n_col),tty(n_col),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if - - t2l(:) = zero - w2l(:) = zero - tx(:) = zero - ty(:) = zero - tty(:) = zero - - ! - ! Need temp copies to handle Y<- betaY + K^-1 X - ! One of the temp copies is not strictly needed when beta==zero - ! - call psb_axpby(one,x,zero,tx,desc_data,info) - call psb_axpby(one,y,zero,ty,desc_data,info) - if(info /=0) goto 9999 - - call psb_baseprc_aply(baseprecv(1),tx,zero,tty,desc_data,trans,work,info) - if(info /=0) goto 9999 - - call psb_spmm(-one,baseprecv(2)%aorig,tty,one,tx,desc_data,info,work=work) - if(info /=0) goto 9999 - - if (ismth /= no_smth_) then - if (baseprecv(2)%iprcparm(glb_smth_) >0) then - call psb_halo(tx,baseprecv(1)%desc_data,info,work=work) - if(info /=0) goto 9999 - else - tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zero - end if - - call psb_csmm(one,baseprecv(2)%av(sm_pr_t_),tx,zero,t2l,info) - if(info /=0) goto 9999 - else - ! - ! Raw aggregation, may take shortcuts - ! - do i=1,desc_data%matrix_data(psb_n_row_) - t2l(baseprecv(2)%mlia(i)) = t2l(baseprecv(2)%mlia(i)) + tx(i) - end do - end if - - - if (baseprecv(2)%iprcparm(coarse_mat_)==mat_repl_) Then - call gsum2d(icontxt,'All',t2l(1:nrg)) - else if (baseprecv(2)%iprcparm(coarse_mat_) /= mat_distr_) Then - write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',& - & baseprecv(2)%iprcparm(coarse_mat_) - endif - - - t6 = mpi_wtime() - w2l=t2l - call psb_baseprc_aply(baseprecv(2),w2l,zero,t2l,baseprecv(2)%desc_data,& - & 'N',work,info) - if(info /=0) goto 9999 - - if (ismth /= no_smth_) then - - if (ismth == smth_omg_) & - & call psb_halo(t2l,baseprecv(2)%desc_data,info,work=work) - call psb_csmm(one,baseprecv(2)%av(sm_pr_),t2l,zero,ty,info) - if(info /=0) goto 9999 - - call psb_axpby(one,ty,one,tty,desc_data,info) - if(info /=0) goto 9999 - - else - - do i=1, desc_data%matrix_data(psb_n_row_) - tty(i) = tty(i) + t2l(baseprecv(2)%mlia(i)) - enddo - - end if - - call psb_axpby(one,x,zero,tx,desc_data,info) - if(info /=0) goto 9999 - - call psb_spmm(-one,baseprecv(2)%aorig,tty,one,tx,desc_data,info,work=work) - if(info /=0) goto 9999 - call psb_baseprc_aply(baseprecv(1),tx,one,tty,desc_data,'N',work,info) - - - call psb_axpby(one,tty,beta,y,desc_data,info) - - deallocate(t2l,w2l,tx,ty,tty) - - - case default - - write(0,*) 'Unknown value for ml_smooth_pos',baseprecv(2)%iprcparm(smth_pos_) - - end select - - case default - write(0,*) me, 'Wrong mltype into PRC_APLY ',& - & baseprecv(2)%iprcparm(ml_type_) - end select - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_errpush(info,name) - call psb_erractionrestore(err_act) - if (err_act.eq.act_abort) then - call psb_error() - return - end if - return - -end subroutine psb_dmlprc_aply - -!!$ -!!$ -!!$ MPcube: Multilevel Parallel Preconditioners Package -!!$ for -!!$ Parallel Sparse BLAS v2.0 -!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata -!!$ Alfredo Buttari University of Rome Tor Vergata -!!$ Daniela Di Serafino II University of Naples -!!$ Pasqua D'Ambra ICAR-CNR -!!$ -!!$ Redistribution and use in source and binary forms, with or without -!!$ modification, are permitted provided that the following conditions -!!$ are met: -!!$ 1. Redistributions of source code must retain the above copyright -!!$ notice, this list of conditions and the following disclaimer. -!!$ 2. Redistributions in binary form must reproduce the above copyright -!!$ notice, this list of conditions, and the following disclaimer in the -!!$ documentation and/or other materials provided with the distribution. -!!$ 3. The name of the MPCUBE group or the names of its contributors may -!!$ not be used to endorse or promote products derived from this -!!$ software without specific written permission. -!!$ -!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR -!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MPCUBE GROUP OR ITS CONTRIBUTORS -!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR -!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF -!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS -!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN -!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -!!$ POSSIBILITY OF SUCH DAMAGE. -!!$ -!!$ - -subroutine psb_dprc_aply1(prec,x,desc_data,info,trans) - - use psb_serial_mod - use psb_descriptor_type - use psb_prec_type - use psb_psblas_mod - use psb_const_mod - use psb_error_mod - implicit none - - type(psb_desc_type),intent(in) :: desc_data - type(psb_dprec_type), intent(in) :: prec - real(kind(0.d0)),intent(inout) :: x(:) - integer, intent(out) :: info - character(len=1), optional :: trans - logical,parameter :: debug=.false., debugprt=.false. - real(kind(1.d0)), parameter :: one=1.d0, zero=0.d0 - - interface - subroutine psb_dprc_aply(prec,x,y,desc_data,info,trans, work) - - use psb_descriptor_type - use psb_prec_type - implicit none - - type(psb_desc_type),intent(in) :: desc_data - type(psb_dprec_type), intent(in) :: prec - real(kind(0.d0)),intent(inout) :: x(:), y(:) - integer, intent(out) :: info - character(len=1), optional :: trans - real(kind(0.d0)), optional, target :: work(:) - end subroutine psb_dprc_aply - end interface - - ! Local variables - character :: trans_ - integer :: icontxt,nprow,npcol,me,mycol,i, isz, err_act, int_err(5) - real(kind(1.d0)), pointer :: WW(:), w1(:) - character(len=20) :: name, ch_err - name='psb_dprec1' - info = 0 - call psb_erractionsave(err_act) - - - icontxt=desc_data%matrix_data(psb_ctxt_) - call blacs_gridinfo(icontxt,nprow,npcol,me,mycol) - if (present(trans)) then - trans_=trans - else - trans_='N' - end if - - allocate(ww(size(x)),w1(size(x)),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if - if (debug) write(0,*) 'Prc_aply1 Size(x) ',size(x), size(ww),size(w1) - call psb_dprc_aply(prec,x,ww,desc_data,info,trans_,work=w1) - if(info /=0) goto 9999 - x(:) = ww(:) - deallocate(ww,W1) - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_errpush(info,name) - call psb_erractionrestore(err_act) - if (err_act.eq.act_abort) then - call psb_error() - return - end if - return -end subroutine psb_dprc_aply1