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

Typographical changes for old/new style comparison operators.

Browse Source 
Taken out debug local parameters.
psblas3-type-indexed
Salvatore Filippone 17 years ago
parent df824ff3be
commit 22af704575
15 changed files with 150 additions and 280 deletions
Show Stats Download Patch File Download Diff File

43
prec/psb_dbjac_aply.f90
Unescape Escape View File

@ -53,20 +53,23 @@ subroutine psb_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
integer :: n_row,n_col
real(kind(1.d0)), pointer :: ww(:), aux(:)
integer :: ictxt,np,me, err_act, int_err(5)
logical,parameter :: debug=.false., debugprt=.false.
character(len=20) :: name, ch_err
integer :: debug_level, debug_unit
character :: trans_
character(len=20) :: name, ch_err
name='psb_bjac_aply'
info = 0
call psb_erractionsave(err_act)
ictxt=psb_cd_get_context(desc_data)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = psb_cd_get_context(desc_data)
call psb_info(ictxt, me, np)
select case(trans)
case('N','n')
case('T','t','C','c')
trans_ = toupper(trans)
select case(trans_)
case('N','T','C')
! Ok
case default
call psb_errpush(40,name)
goto 9999
@ -97,37 +100,35 @@ subroutine psb_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
endif
select case(prec%iprcparm(f_type_))
case(f_ilu_n_,f_ilu_e_)
select case(trans)
case('N','n')
select case(trans_)
case('N')
call psb_spsm(done,prec%av(l_pr_),x,dzero,ww,desc_data,info,&
& trans='N',unit='L',diag=prec%d,choice=psb_none_,work=aux)
& trans=trans_,unit='L',diag=prec%d,choice=psb_none_,work=aux)
if(info /=0) goto 9999
call psb_spsm(alpha,prec%av(u_pr_),ww,beta,y,desc_data,info,&
& trans='N',unit='U',choice=psb_none_, work=aux)
& trans=trans_,unit='U',choice=psb_none_, work=aux)
if(info /=0) goto 9999
case('T','t','C','c')
case('T','C')
call psb_spsm(done,prec%av(u_pr_),x,dzero,ww,desc_data,info,&
& trans=trans,unit='L',diag=prec%d,choice=psb_none_, work=aux)
& trans=trans_,unit='L',diag=prec%d,choice=psb_none_, work=aux)
if(info /=0) goto 9999
call psb_spsm(alpha,prec%av(l_pr_),ww,beta,y,desc_data,info,&
& trans=trans,unit='U',choice=psb_none_,work=aux)
& trans=trans_,unit='U',choice=psb_none_,work=aux)
if(info /=0) goto 9999
end select
case default
write(0,*) 'Unknown factorization type in bjac_aply',prec%iprcparm(f_type_)
info = 4001
call psb_errpush(info,name,a_err='Invalid factorization')
goto 9999
end select
if (debugprt) write(0,*)' Y: ',y(:)
if (n_col <= size(work)) then
if ((4*n_col+n_col) <= size(work)) then
@ -145,7 +146,7 @@ subroutine psb_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
9999 continue
call psb_errpush(info,name,i_err=int_err,a_err=ch_err)
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
if (err_act == psb_act_abort_) then
call psb_error()
return
end if

47
prec/psb_dbjac_bld.f90
Unescape Escape View File

@ -47,13 +47,12 @@ subroutine psb_dbjac_bld(a,desc_a,p,upd,info)
character :: trans, unitd
type(psb_dspmat_type) :: atmp
real(kind(1.d0)) :: t1,t2,t3,t4,t5,t6, t7, t8
logical, parameter :: debugprt=.false., debug=.false., aggr_dump=.false.
integer nztota, err_act, n_row, nrow_a,n_col, nhalo
integer :: ictxt,np,me
character(len=20) :: name, ch_err
if(psb_get_errstatus().ne.0) return
if(psb_get_errstatus() /= 0) return
info=0
name='psb_dbjac_bld'
call psb_erractionsave(err_act)
@ -106,12 +105,9 @@ subroutine psb_dbjac_bld(a,desc_a,p,upd,info)
goto 9999
end if
endif
!!$ call psb_csprt(50+me,a,head='% (A)')
nrow_a = psb_cd_get_local_rows(desc_a)
nztota = psb_sp_get_nnzeros(a)
if (debug) write(0,*)me,': out get_nnzeros',nztota
if (debug) call psb_barrier(ictxt)
n_col = psb_cd_get_local_cols(desc_a)
nhalo = n_col-nrow_a
@ -146,17 +142,6 @@ subroutine psb_dbjac_bld(a,desc_a,p,upd,info)
! This is where we have mo renumbering, thus no need
! for ATMP
if (debugprt) then
open(40+me)
call psb_barrier(ictxt)
call psb_csprt(40+me,a,iv=p%desc_data%loc_to_glob,&
& head='% Local matrix')
close(40+me)
endif
t5= psb_wtime()
if (debug) write(0,*) me,' Going for ilu_fct'
if (debug) call psb_barrier(ictxt)
call psb_ilu_fct(a,p%av(l_pr_),p%av(u_pr_),p%d,info)
if(info/=0) then
info=4010
@ -164,30 +149,6 @@ subroutine psb_dbjac_bld(a,desc_a,p,upd,info)
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if (debug) write(0,*) me,' Done dilu_fct'
if (debugprt) then
!
! Print out the factors on file.
!
open(80+me)
call psb_csprt(80+me,p%av(l_pr_),head='% Local L factor')
write(80+me,*) '% Diagonal: ',p%av(l_pr_)%m
do i=1,p%av(l_pr_)%m
write(80+me,*) i,i,p%d(i)
enddo
call psb_csprt(80+me,p%av(u_pr_),head='% Local U factor')
close(80+me)
end if
! ierr = MPE_Log_event( ifcte, 0, "st SIMPLE" )
t6 = psb_wtime()
!
! write(0,'(i3,1x,a,3(1x,g18.9))') me,'renum/factor time',t3-t2,t6-t5
! if (me==0) write(0,'(a,3(1x,g18.9))') 'renum/factor time',t3-t2,t6-t5
if (psb_sp_getifld(psb_upd_,p%av(u_pr_),info) /= psb_upd_perm_) then
call psb_sp_trim(p%av(u_pr_),info)
@ -199,15 +160,12 @@ subroutine psb_dbjac_bld(a,desc_a,p,upd,info)
case(f_none_)
write(0,*) 'Fact=None in BASEPRC_BLD Bja/ASM??'
info=4010
ch_err='Inconsistent prec f_none_'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
case default
write(0,*) 'Unknown factor type in baseprc_bld bja/asm: ',&
&p%iprcparm(f_type_)
info=4010
ch_err='Unknown f_type_'
call psb_errpush(info,name,a_err=ch_err)
@ -215,13 +173,12 @@ subroutine psb_dbjac_bld(a,desc_a,p,upd,info)
end select
if (debug) write(0,*) me,'End of ilu_bld'
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
if (err_act == psb_act_abort_) then
call psb_error()
return
end if

14
prec/psb_ddiagsc_bld.f90
Unescape Escape View File

@ -46,17 +46,15 @@ subroutine psb_ddiagsc_bld(a,desc_a,p,upd,info)
& me,np,mglob, err_act
integer :: int_err(5)
logical, parameter :: debug=.false.
integer,parameter :: iroot=psb_root_,iout=60,ilout=40
character(len=20) :: name, ch_err
if(psb_get_errstatus().ne.0) return
if(psb_get_errstatus() /= 0) return
info=0
err=0
call psb_erractionsave(err_act)
name = 'psb_diagsc_bld'
if (debug) write(0,*) 'Entering diagsc_bld'
info = 0
int_err(1) = 0
ictxt = psb_cd_get_context(desc_a)
@ -64,10 +62,8 @@ subroutine psb_ddiagsc_bld(a,desc_a,p,upd,info)
n_col = psb_cd_get_local_cols(desc_a)
mglob = psb_cd_get_global_rows(desc_a)
if (debug) write(0,*) 'Preconditioner Blacs_gridinfo'
call psb_info(ictxt, me, np)
if (debug) write(0,*) 'Precond: Diagonal scaling'
! diagonal scaling
call psb_realloc(n_col,p%d,info)
@ -95,7 +91,6 @@ subroutine psb_ddiagsc_bld(a,desc_a,p,upd,info)
goto 9999
end if
if (debug) write(ilout+me,*) 'VDIAG ',n_row
!
! The i-th diagonal entry of the preconditioner is set to one if the
! corresponding entry a_ii of the sparse matrix A is zero; otherwise
@ -107,8 +102,6 @@ subroutine psb_ddiagsc_bld(a,desc_a,p,upd,info)
else
p%d(i) = done/p%d(i)
endif
if (debug) write(ilout+me,*) i,desc_a%loc_to_glob(i), p%d(i)
end do
if (a%pl(1) /= 0) then
@ -124,15 +117,12 @@ subroutine psb_ddiagsc_bld(a,desc_a,p,upd,info)
end if
endif
if (debug) write(*,*) 'Preconditioner DIAG computed OK'
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
if (err_act == psb_act_abort_) then
call psb_error()
return
end if

21
prec/psb_dgprec_aply.f90
Unescape Escape View File

@ -50,9 +50,8 @@ subroutine psb_dgprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
! Local variables
integer :: n_row,int_err(5)
real(kind(1.d0)), pointer :: ww(:)
character ::diagl, diagu
character :: trans_
integer :: ictxt,np,me, err_act
logical,parameter :: debug=.false., debugprt=.false.
character(len=20) :: name, ch_err
name='psb_dgprec_aply'
@ -62,12 +61,11 @@ subroutine psb_dgprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
ictxt=desc_data%matrix_data(psb_ctxt_)
call psb_info(ictxt, me, np)
diagl='U'
diagu='U'
trans_ = toupper(trans)
select case(trans)
case('N','n')
case('T','t','C','c')
select case(trans_)
case('N')
case('T','C')
case default
info=40
int_err(1)=6
@ -108,15 +106,16 @@ subroutine psb_dgprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
case(bjac_)
call psb_bjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
if(info.ne.0) then
if(info /= 0) then
info=4010
ch_err='psb_bjac_aply'
goto 9999
end if
case default
write(0,*) 'Invalid PRE%PREC ',prec%iprcparm(p_type_),':',&
& min_prec_,noprec_,diag_,bjac_
info = 4001
call psb_errpush(info,name,a_err='Invalid prectype')
goto 9999
end select
call psb_erractionrestore(err_act)
@ -125,7 +124,7 @@ subroutine psb_dgprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
9999 continue
call psb_errpush(info,name,i_err=int_err,a_err=ch_err)
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
if (err_act == psb_act_abort_) then
call psb_error()
return
end if

42
prec/psb_dilu_fct.f90
Unescape Escape View File

@ -48,7 +48,6 @@ subroutine psb_dilu_fct(a,l,u,d,info,blck)
integer :: l1,l2,m,err_act
type(psb_dspmat_type), pointer :: blck_
character(len=20) :: name, ch_err
logical, parameter :: debug=.false.
name='psb_dcsrlu'
info = 0
call psb_erractionsave(err_act)
@ -66,7 +65,7 @@ subroutine psb_dilu_fct(a,l,u,d,info,blck)
call psb_nullify_sp(blck_) ! Why do we need this? Who knows....
call psb_sp_all(0,0,blck_,1,info)
if(info.ne.0) then
if(info /= 0) then
info=4010
ch_err='psb_sp_all'
call psb_errpush(info,name,a_err=ch_err)
@ -76,10 +75,9 @@ subroutine psb_dilu_fct(a,l,u,d,info,blck)
blck_%m=0
endif
!!$ write(0,*) 'ilu_fct: ',size(l%ia2),size(u%ia2),a%m,blck_%m
call psb_dilu_fctint(m,a%m,a,blck_%m,blck_,&
& d,l%aspk,l%ia1,l%ia2,u%aspk,u%ia1,u%ia2,l1,l2,info)
if(info.ne.0) then
if(info /= 0) then
info=4010
ch_err='psb_dilu_fctint'
call psb_errpush(info,name,a_err=ch_err)
@ -100,7 +98,7 @@ subroutine psb_dilu_fct(a,l,u,d,info,blck)
blck_ => null()
else
call psb_sp_free(blck_,info)
if(info.ne.0) then
if(info /= 0) then
info=4010
ch_err='psb_sp_free'
call psb_errpush(info,name,a_err=ch_err)
@ -114,7 +112,7 @@ subroutine psb_dilu_fct(a,l,u,d,info,blck)
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
@ -133,37 +131,34 @@ contains
integer :: i,j,k,l,low1,low2,kk,jj,ll, irb, ktrw,err_act
real(kind(1.d0)) :: dia,temp
integer, parameter :: nrb=16
logical,parameter :: debug=.false.
type(psb_dspmat_type) :: trw
integer :: int_err(5)
character(len=20) :: name, ch_err
name='psb_dilu_fctint'
if(psb_get_errstatus().ne.0) return
if(psb_get_errstatus() /= 0) return
info=0
call psb_erractionsave(err_act)
call psb_nullify_sp(trw)
trw%m=0
trw%k=0
if(debug) write(0,*)'LUINT Allocating TRW'
call psb_sp_all(trw,1,info)
if(info.ne.0) then
if(info /= 0) then
info=4010
ch_err='psb_sp_all'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if(debug) write(0,*)'LUINT Done Allocating TRW'
lia2(1) = 1
uia2(1) = 1
l1=0
l2=0
m = ma+mb
if(debug) write(0,*)'In DCSRLU Begin cycle',m,ma,mb
do i = 1, ma
if(debug) write(0,*)'LUINT: Loop index ',i,ma
d(i) = 0.d0
d(i) = dzero
!
! Here we take a fast shortcut if possible, otherwise
@ -192,7 +187,7 @@ contains
if ((mod(i,nrb) == 1).or.(nrb==1)) then
irb = min(ma-i+1,nrb)
call psb_sp_getblk(i,a,trw,info,lrw=i+irb-1)
if(info.ne.0) then
if(info /= 0) then
info=4010
ch_err='psb_sp_getblk'
call psb_errpush(info,name,a_err=ch_err)
@ -295,7 +290,7 @@ contains
call psb_errpush(info,name,i_err=int_err,a_err=ch_err)
goto 9999
else
dia = 1.d0/dia
dia = done/dia
end if
d(i) = dia
! write(6,*)'diag(',i,')=',d(i)
@ -306,19 +301,19 @@ contains
enddo
do i = ma+1, m
d(i) = 0.d0
d(i) = dzero
if (b%fida=='CSR') then
do j = b%ia2(i-ma), b%ia2(i-ma+1) - 1
k = b%ia1(j)
! if (me.eq.2) write(0,*)'ecco k=',k
! if (me == 2) write(0,*)'ecco k=',k
if ((k < i).and.(k >= 1)) then
l1 = l1 + 1
laspk(l1) = b%aspk(j)
lia1(l1) = k
! if(me.eq.2) write(0,*)'scrivo l'
! if(me == 2) write(0,*)'scrivo l'
else if (k == i) then
d(i) = b%aspk(j)
else if ((k > i).and.(k <= m)) then
@ -334,7 +329,7 @@ contains
if ((mod((i-ma),nrb) == 1).or.(nrb==1)) then
irb = min(m-i+1,nrb)
call psb_sp_getblk(i-ma,b,trw,info,lrw=i-ma+irb-1)
if(info.ne.0) then
if(info /= 0) then
info=4010
ch_err='psb_sp_getblk'
call psb_errpush(info,name,a_err=ch_err)
@ -435,7 +430,7 @@ contains
call psb_errpush(info,name,i_err=int_err,a_err=ch_err)
goto 9999
else
dia = 1.d0/dia
dia = done/dia
end if
d(i) = dia
! Scale row i of upper triangle
@ -445,20 +440,19 @@ contains
enddo
call psb_sp_free(trw,info)
if(info.ne.0) then
if(info /= 0) then
info=4010
ch_err='psb_sp_free'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if(debug) write(0,*)'Leaving ilu_fct'
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
if (err_act == psb_act_abort_) then
call psb_error()
return
end if

12
prec/psb_dprc_aply.f90
Unescape Escape View File

@ -46,7 +46,6 @@ subroutine psb_dprc_aply(prec,x,y,desc_data,info,trans, work)
character :: trans_
real(kind(1.d0)), pointer :: work_(:)
integer :: ictxt,np,me,err_act
logical,parameter :: debug=.false., debugprt=.false.
character(len=20) :: name
name='psb_prc_aply'
@ -74,7 +73,10 @@ subroutine psb_dprc_aply(prec,x,y,desc_data,info,trans, work)
end if
call psb_gprec_aply(done,prec,x,dzero,y,desc_data,trans_,work_,info)
! If the original distribution has an overlap we should fix that.
call psb_ovrl(y,desc_data,info,update=psb_avg_)
if (present(work)) then
else
deallocate(work_)
@ -85,7 +87,7 @@ subroutine psb_dprc_aply(prec,x,y,desc_data,info,trans, work)
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
@ -141,7 +143,6 @@ subroutine psb_dprc_aply1(prec,x,desc_data,info,trans)
real(kind(0.d0)),intent(inout) :: x(:)
integer, intent(out) :: info
character(len=1), optional :: trans
logical,parameter :: debug=.false., debugprt=.false.
! Local variables
character :: trans_
@ -156,7 +157,7 @@ subroutine psb_dprc_aply1(prec,x,desc_data,info,trans)
ictxt=desc_data%matrix_data(psb_ctxt_)
call psb_info(ictxt, me, np)
if (present(trans)) then
trans_=trans
trans_=toupper(trans)
else
trans_='N'
end if
@ -166,7 +167,6 @@ subroutine psb_dprc_aply1(prec,x,desc_data,info,trans)
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(:)
@ -178,7 +178,7 @@ subroutine psb_dprc_aply1(prec,x,desc_data,info,trans)
9999 continue
call psb_errpush(info,name)
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
if (err_act == psb_act_abort_) then
call psb_error()
return
end if

30
prec/psb_dprecbld.f90
Unescape Escape View File

@ -44,35 +44,32 @@ subroutine psb_dprecbld(a,desc_a,p,info,upd)
Integer :: err, n_row, n_col,ictxt,&
& me,np,mglob, err_act
integer :: int_err(5)
character :: iupd
character :: upd_
logical, parameter :: debug=.false.
integer,parameter :: iroot=psb_root_,iout=60,ilout=40
character(len=20) :: name, ch_err
if(psb_get_errstatus().ne.0) return
if(psb_get_errstatus() /= 0) return
info=0
err=0
call psb_erractionsave(err_act)
name = 'psb_precbld'
if (debug) write(0,*) 'Entering precbld',P%prec,desc_a%matrix_data(:)
info = 0
int_err(1) = 0
ictxt = psb_cd_get_context(desc_a)
if (debug) write(0,*) 'Preconditioner psb_info'
call psb_info(ictxt, me, np)
if (present(upd)) then
if (debug) write(0,*) 'UPD ', upd
if ((upd.eq.'F').or.(upd.eq.'T')) then
iupd=upd
else
iupd='F'
endif
upd_ = toupper(upd)
else
iupd='F'
upd_='F'
endif
if ((upd_ == 'F').or.(upd_ == 'T')) then
! ok
else
upd_='F'
endif
n_row = psb_cd_get_local_rows(desc_a)
n_col = psb_cd_get_local_cols(desc_a)
@ -101,8 +98,7 @@ subroutine psb_dprecbld(a,desc_a,p,info,upd)
case (diag_)
call psb_diagsc_bld(a,desc_a,p,iupd,info)
if(debug) write(0,*)me,': out of psb_diagsc_bld'
call psb_diagsc_bld(a,desc_a,p,upd_,info)
if(info /= 0) then
info=4010
ch_err='psb_diagsc_bld'
@ -115,9 +111,7 @@ subroutine psb_dprecbld(a,desc_a,p,info,upd)
call psb_check_def(p%iprcparm(f_type_),'fact',&
& f_ilu_n_,is_legal_ml_fact)
if (debug) write(0,*)me, ': Calling PSB_BJAC_BLD'
if (debug) call psb_barrier(ictxt)
call psb_bjac_bld(a,desc_a,p,iupd,info)
call psb_bjac_bld(a,desc_a,p,upd_,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_bjac_bld')
@ -137,7 +131,7 @@ subroutine psb_dprecbld(a,desc_a,p,info,upd)
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
if (err_act == psb_act_abort_) then
call psb_error()
return
end if

8
prec/psb_prec_type.f90
Unescape Escape View File

@ -226,7 +226,7 @@ contains
integer, intent(out) :: info
integer :: me, err_act,i
character(len=20) :: name
if(psb_get_errstatus().ne.0) return
if(psb_get_errstatus() /= 0) return
info=0
name = 'psb_precfree'
call psb_erractionsave(err_act)
@ -277,7 +277,7 @@ contains
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
@ -300,7 +300,7 @@ contains
integer, intent(out) :: info
integer :: err_act,i
character(len=20) :: name
if(psb_get_errstatus().ne.0) return
if(psb_get_errstatus() /= 0) return
info=0
name = 'psb_precfree'
call psb_erractionsave(err_act)
@ -345,7 +345,7 @@ contains
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
if (err_act == psb_act_abort_) then
call psb_error()
return
end if

53
prec/psb_zbjac_aply.f90
Unescape Escape View File

@ -53,20 +53,23 @@ subroutine psb_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
integer :: n_row,n_col
complex(kind(1.d0)), pointer :: ww(:), aux(:)
integer :: ictxt,np,me, err_act, int_err(5)
logical,parameter :: debug=.false., debugprt=.false.
character(len=20) :: name, ch_err
integer :: debug_level, debug_unit
character :: trans_
character(len=20) :: name, ch_err
name='psb_bjac_aply'
info = 0
call psb_erractionsave(err_act)
ictxt=psb_cd_get_context(desc_data)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = psb_cd_get_context(desc_data)
call psb_info(ictxt, me, np)
select case(trans)
case('N','n')
case('T','t','C','c')
trans_ = toupper(trans)
select case(trans_)
case('N','T','C')
! Ok
case default
call psb_errpush(40,name)
goto 9999
@ -97,37 +100,36 @@ subroutine psb_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
endif
select case(prec%iprcparm(f_type_))
case(f_ilu_n_,f_ilu_e_)
select case(trans)
case('N','n')
select case(trans_)
case('N')
call psb_spsm(zone,prec%av(l_pr_),x,zzero,ww,desc_data,info,&
& trans='N',unit='L',diag=prec%d,choice=psb_none_,work=aux)
& trans=trans_,unit='L',diag=prec%d,choice=psb_none_,work=aux)
if(info /=0) goto 9999
call psb_spsm(alpha,prec%av(u_pr_),ww,beta,y,desc_data,info,&
& trans='N',unit='U',choice=psb_none_, work=aux)
& trans=trans_,unit='U',choice=psb_none_, work=aux)
if(info /=0) goto 9999
case('T','t','C','c')
case('T','C')
call psb_spsm(zone,prec%av(u_pr_),x,zzero,ww,desc_data,info,&
& trans=trans,unit='L',diag=prec%d,choice=psb_none_, work=aux)
& trans=trans_,unit='L',diag=prec%d,choice=psb_none_, work=aux)
if(info /=0) goto 9999
call psb_spsm(alpha,prec%av(l_pr_),ww,beta,y,desc_data,info,&
& trans=trans,unit='U',choice=psb_none_,work=aux)
& trans=trans_,unit='U',choice=psb_none_,work=aux)
if(info /=0) goto 9999
end select
case default
write(0,*) 'Unknown factorization type in bjac_aply',prec%iprcparm(f_type_)
info = 4001
call psb_errpush(info,name,a_err='Invalid factorization')
goto 9999
end select
if (debugprt) write(0,*)' Y: ',y(:)
if (n_col <= size(work)) then
if ((4*n_col+n_col) <= size(work)) then
else
@ -144,11 +146,10 @@ subroutine psb_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
9999 continue
call psb_errpush(info,name,i_err=int_err,a_err=ch_err)
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_zbjac_aply

47
prec/psb_zbjac_bld.f90
Unescape Escape View File

@ -47,13 +47,12 @@ subroutine psb_zbjac_bld(a,desc_a,p,upd,info)
character :: trans, unitd
type(psb_zspmat_type) :: atmp
real(kind(1.d0)) :: t1,t2,t3,t4,t5,t6, t7, t8
logical, parameter :: debugprt=.false., debug=.false., aggr_dump=.false.
integer nztota, err_act, n_row, nrow_a,n_col, nhalo
integer :: ictxt,np,me
character(len=20) :: name, ch_err
if(psb_get_errstatus().ne.0) return
if(psb_get_errstatus() /= 0) return
info=0
name='psb_zbjac_bld'
call psb_erractionsave(err_act)
@ -106,12 +105,9 @@ subroutine psb_zbjac_bld(a,desc_a,p,upd,info)
goto 9999
end if
endif
!!$ call psb_csprt(50+me,a,head='% (A)')
nrow_a = psb_cd_get_local_rows(desc_a)
nztota = psb_sp_get_nnzeros(a)
if (debug) write(0,*)me,': out get_nnzeros',nztota
if (debug) call psb_barrier(ictxt)
n_col = psb_cd_get_local_cols(desc_a)
nhalo = n_col-nrow_a
@ -146,17 +142,6 @@ subroutine psb_zbjac_bld(a,desc_a,p,upd,info)
! This is where we have mo renumbering, thus no need
! for ATMP
if (debugprt) then
open(40+me)
call psb_barrier(ictxt)
call psb_csprt(40+me,a,iv=p%desc_data%loc_to_glob,&
& head='% Local matrix')
close(40+me)
endif
t5= psb_wtime()
if (debug) write(0,*) me,' Going for ilu_fct'
if (debug) call psb_barrier(ictxt)
call psb_ilu_fct(a,p%av(l_pr_),p%av(u_pr_),p%d,info)
if(info/=0) then
info=4010
@ -164,30 +149,6 @@ subroutine psb_zbjac_bld(a,desc_a,p,upd,info)
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if (debug) write(0,*) me,' Done dilu_fct'
if (debugprt) then
!
! Print out the factors on file.
!
open(80+me)
call psb_csprt(80+me,p%av(l_pr_),head='% Local L factor')
write(80+me,*) '% Diagonal: ',p%av(l_pr_)%m
do i=1,p%av(l_pr_)%m
write(80+me,*) i,i,p%d(i)
enddo
call psb_csprt(80+me,p%av(u_pr_),head='% Local U factor')
close(80+me)
end if
! ierr = MPE_Log_event( ifcte, 0, "st SIMPLE" )
t6 = psb_wtime()
!
! write(0,'(i3,1x,a,3(1x,g18.9))') me,'renum/factor time',t3-t2,t6-t5
! if (me==0) write(0,'(a,3(1x,g18.9))') 'renum/factor time',t3-t2,t6-t5
if (psb_sp_getifld(psb_upd_,p%av(u_pr_),info) /= psb_upd_perm_) then
call psb_sp_trim(p%av(u_pr_),info)
@ -199,15 +160,12 @@ subroutine psb_zbjac_bld(a,desc_a,p,upd,info)
case(f_none_)
write(0,*) 'Fact=None in BASEPRC_BLD Bja/ASM??'
info=4010
ch_err='Inconsistent prec f_none_'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
case default
write(0,*) 'Unknown factor type in baseprc_bld bja/asm: ',&
&p%iprcparm(f_type_)
info=4010
ch_err='Unknown f_type_'
call psb_errpush(info,name,a_err=ch_err)
@ -215,13 +173,12 @@ subroutine psb_zbjac_bld(a,desc_a,p,upd,info)
end select
if (debug) write(0,*) me,'End of ilu_bld'
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
if (err_act == psb_act_abort_) then
call psb_error()
return
end if

14
prec/psb_zdiagsc_bld.f90
Unescape Escape View File

@ -46,17 +46,15 @@ subroutine psb_zdiagsc_bld(a,desc_a,p,upd,info)
& me,np,mglob,err_act
integer :: int_err(5)
logical, parameter :: debug=.false.
integer,parameter :: iroot=psb_root_,iout=60,ilout=40
character(len=20) :: name, ch_err
if(psb_get_errstatus().ne.0) return
if(psb_get_errstatus() /= 0) return
info=0
err=0
call psb_erractionsave(err_act)
name = 'psb_diagsc_bld'
if (debug) write(0,*) 'Entering diagsc_bld'
info = 0
int_err(1) = 0
ictxt = psb_cd_get_context(desc_a)
@ -64,10 +62,8 @@ subroutine psb_zdiagsc_bld(a,desc_a,p,upd,info)
n_col = psb_cd_get_local_cols(desc_a)
mglob = psb_cd_get_global_rows(desc_a)
if (debug) write(0,*) 'Preconditioner Blacs_gridinfo'
call psb_info(ictxt, me, np)
if (debug) write(0,*) 'Precond: Diagonal scaling'
! diagonal scaling
call psb_realloc(n_col,p%d,info)
@ -95,7 +91,6 @@ subroutine psb_zdiagsc_bld(a,desc_a,p,upd,info)
goto 9999
end if
if (debug) write(ilout+me,*) 'VDIAG ',n_row
!
! The i-th diagonal entry of the preconditioner is set to one if the
! corresponding entry a_ii of the sparse matrix A is zero; otherwise
@ -107,8 +102,6 @@ subroutine psb_zdiagsc_bld(a,desc_a,p,upd,info)
else
p%d(i) = zone/p%d(i)
endif
if (debug) write(ilout+me,*) i,desc_a%loc_to_glob(i), p%d(i)
end do
if (a%pl(1) /= 0) then
@ -124,15 +117,12 @@ subroutine psb_zdiagsc_bld(a,desc_a,p,upd,info)
end if
endif
if (debug) write(*,*) 'Preconditioner DIAG computed OK'
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
if (err_act == psb_act_abort_) then
call psb_error()
return
end if

21
prec/psb_zgprec_aply.f90
Unescape Escape View File

@ -50,9 +50,8 @@ subroutine psb_zgprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
! Local variables
integer :: n_row,int_err(5)
complex(kind(1.d0)), pointer :: ww(:)
character ::diagl, diagu
character :: trans_
integer :: ictxt,np,me, err_act
logical,parameter :: debug=.false., debugprt=.false.
character(len=20) :: name, ch_err
@ -63,12 +62,11 @@ subroutine psb_zgprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
ictxt=desc_data%matrix_data(psb_ctxt_)
call psb_info(ictxt, me, np)
diagl='U'
diagu='U'
trans_ = toupper(trans)
select case(trans)
case('N','n')
case('T','t','C','c')
select case(trans_)
case('N')
case('T','C')
case default
info=40
int_err(1)=6
@ -109,15 +107,16 @@ subroutine psb_zgprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
case(bjac_)
call psb_bjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
if(info.ne.0) then
if(info /= 0) then
info=4010
ch_err='psb_bjac_aply'
goto 9999
end if
case default
write(0,*) 'Invalid PRE%PREC ',prec%iprcparm(p_type_),':',&
& min_prec_,noprec_,diag_,bjac_
info = 4001
call psb_errpush(info,name,a_err='Invalid prectype')
goto 9999
end select
call psb_erractionrestore(err_act)
@ -126,7 +125,7 @@ subroutine psb_zgprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
9999 continue
call psb_errpush(info,name,i_err=int_err,a_err=ch_err)
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
if (err_act == psb_act_abort_) then
call psb_error()
return
end if

32
prec/psb_zilu_fct.f90
Unescape Escape View File

@ -65,7 +65,7 @@ subroutine psb_zilu_fct(a,l,u,d,info,blck)
call psb_nullify_sp(blck_) ! Why do we need this? Who knows....
call psb_sp_all(0,0,blck_,1,info)
if(info.ne.0) then
if(info /= 0) then
info=4010
ch_err='psb_sp_all'
call psb_errpush(info,name,a_err=ch_err)
@ -77,7 +77,7 @@ subroutine psb_zilu_fct(a,l,u,d,info,blck)
call psb_zilu_fctint(m,a%m,a,blck_%m,blck_,&
& d,l%aspk,l%ia1,l%ia2,u%aspk,u%ia1,u%ia2,l1,l2,info)
if(info.ne.0) then
if(info /= 0) then
info=4010
ch_err='psb_zilu_fctint'
call psb_errpush(info,name,a_err=ch_err)
@ -98,7 +98,7 @@ subroutine psb_zilu_fct(a,l,u,d,info,blck)
blck_ => null()
else
call psb_sp_free(blck_,info)
if(info.ne.0) then
if(info /= 0) then
info=4010
ch_err='psb_sp_free'
call psb_errpush(info,name,a_err=ch_err)
@ -112,7 +112,7 @@ subroutine psb_zilu_fct(a,l,u,d,info,blck)
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
@ -131,36 +131,33 @@ contains
integer :: i,j,k,l,low1,low2,kk,jj,ll, irb, ktrw,err_act
complex(kind(1.d0)) :: dia,temp
integer, parameter :: nrb=16
logical,parameter :: debug=.false.
type(psb_zspmat_type) :: trw
integer :: int_err(5)
character(len=20) :: name, ch_err
name='psb_zilu_fctint'
if(psb_get_errstatus().ne.0) return
if(psb_get_errstatus() /= 0) return
info=0
call psb_erractionsave(err_act)
call psb_nullify_sp(trw)
trw%m=0
trw%k=0
if(debug) write(0,*)'LUINT Allocating TRW'
call psb_sp_all(trw,1,info)
if(info.ne.0) then
if(info /= 0) then
info=4010
ch_err='psb_sp_all'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if(debug) write(0,*)'LUINT Done Allocating TRW'
lia2(1) = 1
uia2(1) = 1
l1=0
l2=0
m = ma+mb
if(debug) write(0,*)'In DCSRLU Begin cycle',m,ma,mb
do i = 1, ma
if(debug) write(0,*)'LUINT: Loop index ',i,ma
d(i) = zzero
!
@ -190,7 +187,7 @@ contains
if ((mod(i,nrb) == 1).or.(nrb==1)) then
irb = min(ma-i+1,nrb)
call psb_sp_getblk(i,a,trw,info,lrw=i+irb-1)
if(info.ne.0) then
if(info /= 0) then
info=4010
ch_err='psb_sp_getblk'
call psb_errpush(info,name,a_err=ch_err)
@ -311,12 +308,12 @@ contains
do j = b%ia2(i-ma), b%ia2(i-ma+1) - 1
k = b%ia1(j)
! if (me.eq.2) write(0,*)'ecco k=',k
! if (me == 2) write(0,*)'ecco k=',k
if ((k < i).and.(k >= 1)) then
l1 = l1 + 1
laspk(l1) = b%aspk(j)
lia1(l1) = k
! if(me.eq.2) write(0,*)'scrivo l'
! if(me == 2) write(0,*)'scrivo l'
else if (k == i) then
d(i) = b%aspk(j)
else if ((k > i).and.(k <= m)) then
@ -332,7 +329,7 @@ contains
if ((mod((i-ma),nrb) == 1).or.(nrb==1)) then
irb = min(m-i+1,nrb)
call psb_sp_getblk(i-ma,b,trw,info,lrw=i-ma+irb-1)
if(info.ne.0) then
if(info /= 0) then
info=4010
ch_err='psb_sp_getblk'
call psb_errpush(info,name,a_err=ch_err)
@ -443,20 +440,19 @@ contains
enddo
call psb_sp_free(trw,info)
if(info.ne.0) then
if(info /= 0) then
info=4010
ch_err='psb_sp_free'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if(debug) write(0,*)'Leaving ilu_fct'
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
if (err_act == psb_act_abort_) then
call psb_error()
return
end if

15
prec/psb_zprc_aply.f90
Unescape Escape View File

@ -47,7 +47,6 @@ subroutine psb_zprc_aply(prec,x,y,desc_data,info,trans, work)
character :: trans_
complex(kind(1.d0)), pointer :: work_(:)
integer :: ictxt,np,me,err_act
logical,parameter :: debug=.false., debugprt=.false.
character(len=20) :: name
name='psb_zprec_aply'
@ -73,8 +72,10 @@ subroutine psb_zprc_aply(prec,x,y,desc_data,info,trans, work)
end if
end if
call psb_gprec_aply(zone,prec,x,zzero,y,desc_data,trans_, work_,info)
call psb_gprec_aply(zone,prec,x,zzero,y,desc_data,trans_,work_,info)
! If the original distribution has an overlap we should fix that.
call psb_ovrl(y,desc_data,info,update=psb_avg_)
if (present(work)) then
@ -87,7 +88,7 @@ subroutine psb_zprc_aply(prec,x,y,desc_data,info,trans, work)
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
@ -143,7 +144,6 @@ subroutine psb_zprc_aply1(prec,x,desc_data,info,trans)
complex(kind(0.d0)),intent(inout) :: x(:)
integer, intent(out) :: info
character(len=1), optional :: trans
logical,parameter :: debug=.false., debugprt=.false.
! Local variables
character :: trans_
@ -158,7 +158,7 @@ subroutine psb_zprc_aply1(prec,x,desc_data,info,trans)
ictxt=desc_data%matrix_data(psb_ctxt_)
call psb_info(ictxt, me, np)
if (present(trans)) then
trans_=trans
trans_=toupper(trans)
else
trans_='N'
end if
@ -168,7 +168,6 @@ subroutine psb_zprc_aply1(prec,x,desc_data,info,trans)
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_zprc_aply(prec,x,ww,desc_data,info,trans_,work=w1)
if(info /=0) goto 9999
x(:) = ww(:)
@ -180,7 +179,7 @@ subroutine psb_zprc_aply1(prec,x,desc_data,info,trans)
9999 continue
call psb_errpush(info,name)
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
if (err_act == psb_act_abort_) then
call psb_error()
return
end if

31
prec/psb_zprecbld.f90
Unescape Escape View File

@ -42,39 +42,35 @@ subroutine psb_zprecbld(a,desc_a,p,info,upd)
! Local scalars
Integer :: err, n_row, n_col,ictxt,&
& me,np,mglob, err_act
integer :: int_err(5)
character :: iupd
character :: upd_
logical, parameter :: debug=.false.
integer,parameter :: iroot=psb_root_,iout=60,ilout=40
character(len=20) :: name, ch_err
if(psb_get_errstatus().ne.0) return
if(psb_get_errstatus() /= 0) return
info=0
err=0
call psb_erractionsave(err_act)
name = 'psb_precbld'
if (debug) write(0,*) 'Entering precbld',P%prec,desc_a%matrix_data(:)
info = 0
int_err(1) = 0
ictxt = psb_cd_get_context(desc_a)
if (debug) write(0,*) 'Preconditioner psb_info'
call psb_info(ictxt, me, np)
if (present(upd)) then
if (debug) write(0,*) 'UPD ', upd
if ((upd.eq.'F').or.(upd.eq.'T')) then
iupd=upd
else
iupd='F'
endif
upd_ = toupper(upd)
else
upd_='F'
endif
if ((upd_ == 'F').or.(upd_ == 'T')) then
! ok
else
iupd='F'
upd_='F'
endif
n_row = psb_cd_get_local_rows(desc_a)
n_col = psb_cd_get_local_cols(desc_a)
@ -103,8 +99,7 @@ subroutine psb_zprecbld(a,desc_a,p,info,upd)
case (diag_)
call psb_diagsc_bld(a,desc_a,p,iupd,info)
if(debug) write(0,*)me,': out of psb_diagsc_bld'
call psb_diagsc_bld(a,desc_a,p,upd_,info)
if(info /= 0) then
info=4010
ch_err='psb_diagsc_bld'
@ -117,9 +112,7 @@ subroutine psb_zprecbld(a,desc_a,p,info,upd)
call psb_check_def(p%iprcparm(f_type_),'fact',&
& f_ilu_n_,is_legal_ml_fact)
if (debug) write(0,*)me, ': Calling PSB_BJAC_BLD'
if (debug) call psb_barrier(ictxt)
call psb_bjac_bld(a,desc_a,p,iupd,info)
call psb_bjac_bld(a,desc_a,p,upd_,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_bjac_bld')
@ -139,7 +132,7 @@ subroutine psb_zprecbld(a,desc_a,p,info,upd)
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
if (err_act == psb_act_abort_) then
call psb_error()
return
end if

Write Preview
Loading…
Cancel
Save