diff --git a/prec/psb_dbjac_aply.f90 b/prec/psb_dbjac_aply.f90 index a7c63b68..dbb70af3 100644 --- a/prec/psb_dbjac_aply.f90 +++ b/prec/psb_dbjac_aply.f90 @@ -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 diff --git a/prec/psb_dbjac_bld.f90 b/prec/psb_dbjac_bld.f90 index db67dabe..aae49e63 100644 --- a/prec/psb_dbjac_bld.f90 +++ b/prec/psb_dbjac_bld.f90 @@ -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 diff --git a/prec/psb_ddiagsc_bld.f90 b/prec/psb_ddiagsc_bld.f90 index a5b8337d..7dac53f1 100644 --- a/prec/psb_ddiagsc_bld.f90 +++ b/prec/psb_ddiagsc_bld.f90 @@ -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 diff --git a/prec/psb_dgprec_aply.f90 b/prec/psb_dgprec_aply.f90 index 1d3985fb..30205867 100644 --- a/prec/psb_dgprec_aply.f90 +++ b/prec/psb_dgprec_aply.f90 @@ -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 diff --git a/prec/psb_dilu_fct.f90 b/prec/psb_dilu_fct.f90 index 38a002ea..3b3cb5ec 100644 --- a/prec/psb_dilu_fct.f90 +++ b/prec/psb_dilu_fct.f90 @@ -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 diff --git a/prec/psb_dprc_aply.f90 b/prec/psb_dprc_aply.f90 index e813ec06..1e6f2c1a 100644 --- a/prec/psb_dprc_aply.f90 +++ b/prec/psb_dprc_aply.f90 @@ -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 diff --git a/prec/psb_dprecbld.f90 b/prec/psb_dprecbld.f90 index 304924ea..9ba26257 100644 --- a/prec/psb_dprecbld.f90 +++ b/prec/psb_dprecbld.f90 @@ -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 diff --git a/prec/psb_prec_type.f90 b/prec/psb_prec_type.f90 index 79308dd3..5c95601a 100644 --- a/prec/psb_prec_type.f90 +++ b/prec/psb_prec_type.f90 @@ -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 diff --git a/prec/psb_zbjac_aply.f90 b/prec/psb_zbjac_aply.f90 index 917febe5..c6cb1dd5 100644 --- a/prec/psb_zbjac_aply.f90 +++ b/prec/psb_zbjac_aply.f90 @@ -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 - diff --git a/prec/psb_zbjac_bld.f90 b/prec/psb_zbjac_bld.f90 index fd49f714..7246839f 100644 --- a/prec/psb_zbjac_bld.f90 +++ b/prec/psb_zbjac_bld.f90 @@ -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 diff --git a/prec/psb_zdiagsc_bld.f90 b/prec/psb_zdiagsc_bld.f90 index 71e0d46e..cdfa190d 100644 --- a/prec/psb_zdiagsc_bld.f90 +++ b/prec/psb_zdiagsc_bld.f90 @@ -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 diff --git a/prec/psb_zgprec_aply.f90 b/prec/psb_zgprec_aply.f90 index d02df12d..dc2c4a01 100644 --- a/prec/psb_zgprec_aply.f90 +++ b/prec/psb_zgprec_aply.f90 @@ -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 diff --git a/prec/psb_zilu_fct.f90 b/prec/psb_zilu_fct.f90 index c28901f8..d50206de 100644 --- a/prec/psb_zilu_fct.f90 +++ b/prec/psb_zilu_fct.f90 @@ -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 diff --git a/prec/psb_zprc_aply.f90 b/prec/psb_zprc_aply.f90 index 3e4734ca..4c862149 100644 --- a/prec/psb_zprc_aply.f90 +++ b/prec/psb_zprc_aply.f90 @@ -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 diff --git a/prec/psb_zprecbld.f90 b/prec/psb_zprecbld.f90 index 36a20ab6..06313e39 100644 --- a/prec/psb_zprecbld.f90 +++ b/prec/psb_zprecbld.f90 @@ -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