fdurastante
/
psblas3
mirror of https://github.com/sfilippone/psblas3.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Split the prc_aply section.

Browse Source
psblas3-type-indexed
Salvatore Filippone 19 years ago
parent 580605c087
commit 5f6c77d4a3
6 changed files with 1389 additions and 1386 deletions
Show Stats Download Patch File Download Diff File

7
src/prec/Makefile
Unescape Escape View File

@ -6,8 +6,11 @@ LIBDIR=../../lib/
MPFOBJS=psb_dilu_bld.o psb_dbldaggrmat.o MPFOBJS=psb_dilu_bld.o psb_dbldaggrmat.o
F90OBJS=psb_dasmatbld.o psb_dslu_bld.o psb_dumf_bld.o psb_dilu_fct.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_dmlprc_bld.o psb_dsp_renum.o\
psb_dprec.o psb_dprecbld.o gps.o psb_dprecfree.o psb_dprecset.o \ psb_dprecbld.o gps.o psb_dprecfree.o psb_dprecset.o \
psb_dbaseprc_bld.o psb_dgenaggrmap.o $(MPFOBJS) 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 COBJS=psb_slu_impl.o psb_umf_impl.o
INCDIRS=-I. -I.. -I$(LIBDIR) INCDIRS=-I. -I.. -I$(LIBDIR)

296
src/prec/psb_dbaseprc_aply.f90
Unescape Escape View File

@ -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

298
src/prec/psb_dbjac_aply.f90
Unescape Escape View File

@ -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

529
src/prec/psb_dmlprc_aply.f90
Unescape Escape View File

@ -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

261
src/prec/psb_dprc_aply.f90
Unescape Escape View File

@ -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

1384
src/prec/psb_dprec.f90
View File

File diff suppressed because it is too large Load Diff
Write Preview
Loading…
Cancel
Save