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

Fixed constant names with MLD_ prefix.

Browse Source
stopcriterion
Salvatore Filippone 18 years ago
parent 9f34a33b41
commit 47c27d7e64
44 changed files with 997 additions and 998 deletions
Show Stats Download Patch File Download Diff File

3
Makefile
Unescape Escape View File

@ -10,9 +10,8 @@ kryl:
veryclean:
(cd mlprec; make veryclean)
(cd krylov; make veryclean)
/bin/rm -f $(OBJS) $(LOCAL_MODS)
(cd lib; /bin/rm -f *.a *$(.mod))
clean:
(cd mlprec; make clean)
(cd krylov; make clean)
/bin/rm -f $(OBJS) $(LOCAL_MODS)

6
mlprec/mld_daggrmap_bld.f90
Unescape Escape View File

@ -72,7 +72,7 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
ncol = psb_cd_get_local_cols(desc_a)
select case (aggr_type)
case (dec_aggr_,sym_dec_aggr_)
case (mld_dec_aggr_,sym_mld_dec_aggr_)
nr = a%m
@ -87,7 +87,7 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
do i=1, nr
ilaggr(i) = -(nr+1)
end do
if (aggr_type == dec_aggr_) then
if (aggr_type == mld_dec_aggr_) then
apnt => a
else
call psb_sp_clip(a,atmp,info,imax=nr,jmax=nr,&
@ -305,7 +305,7 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
nlaggr(me+1) = naggr
call psb_sum(ictxt,nlaggr(1:np))
if (aggr_type == sym_dec_aggr_) then
if (aggr_type == sym_mld_dec_aggr_) then
call psb_sp_free(atmp,info)
end if

6
mlprec/mld_daggrmat_asb.f90
Unescape Escape View File

@ -61,8 +61,8 @@ subroutine mld_daggrmat_asb(a,desc_a,ac,desc_ac,p,info)
call psb_info(ictxt, me, np)
select case (p%iprcparm(aggr_kind_))
case (no_smooth_)
select case (p%iprcparm(mld_aggr_kind_))
case (mld_no_smooth_)
call mld_aggrmat_raw_asb(a,desc_a,ac,desc_ac,p,info)
@ -72,7 +72,7 @@ subroutine mld_daggrmat_asb(a,desc_a,ac,desc_ac,p,info)
end if
if (aggr_dump) call psb_csprt(90+me,ac,head='% Raw aggregate.')
case(smooth_prol_,biz_prol_)
case(mld_smooth_prol_,mld_biz_prol_)
if (aggr_dump) call psb_csprt(70+me,a,head='% Input matrix')
call mld_aggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)

8
mlprec/mld_daggrmat_raw_asb.F90
Unescape Escape View File

@ -87,7 +87,7 @@ subroutine mld_daggrmat_raw_asb(a,desc_a,ac,desc_ac,p,info)
naggrm1=sum(p%nlaggr(1:me))
if (p%iprcparm(coarse_mat_) == repl_mat_) then
if (p%iprcparm(mld_coarse_mat_) == mld_repl_mat_) then
do i=1, nrow
p%mlia(i) = p%mlia(i) + naggrm1
end do
@ -120,7 +120,7 @@ subroutine mld_daggrmat_raw_asb(a,desc_a,ac,desc_ac,p,info)
! This is to minimize data exchange
call psb_spcnv(b,info,afmt='coo',dupl=psb_dupl_add_)
if (p%iprcparm(coarse_mat_) == repl_mat_) then
if (p%iprcparm(mld_coarse_mat_) == mld_repl_mat_) then
call psb_cdrep(ntaggr,ictxt,desc_ac,info)
if(info /= 0) then
@ -167,7 +167,7 @@ subroutine mld_daggrmat_raw_asb(a,desc_a,ac,desc_ac,p,info)
goto 9999
end if
else if (p%iprcparm(coarse_mat_) == distr_mat_) then
else if (p%iprcparm(mld_coarse_mat_) == mld_distr_mat_) then
call psb_cdall(ictxt,desc_ac,info,nl=naggr)
if(info /= 0) then
@ -193,7 +193,7 @@ subroutine mld_daggrmat_raw_asb(a,desc_a,ac,desc_ac,p,info)
else
write(0,*) 'Unknown p%iprcparm(coarse_mat) in aggregate_sp',p%iprcparm(coarse_mat_)
write(0,*) 'Unknown p%iprcparm(coarse_mat) in aggregate_sp',p%iprcparm(mld_coarse_mat_)
end if
deallocate(nzbr,idisp)

46
mlprec/mld_daggrmat_smth_asb.F90
Unescape Escape View File

@ -84,7 +84,7 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
call psb_nullify_sp(am3)
call psb_nullify_sp(am4)
am2 => p%av(sm_pr_t_)
am2 => p%av(mld_sm_pr_t_)
am1 => p%av(sm_pr_)
call psb_nullify_sp(am1)
call psb_nullify_sp(am2)
@ -108,9 +108,9 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
naggrm1 = sum(p%nlaggr(1:me))
naggrp1 = sum(p%nlaggr(1:me+1))
ml_global_nmb = ( (p%iprcparm(aggr_kind_) == smooth_prol_).or.&
& ( (p%iprcparm(aggr_kind_) == biz_prol_).and.&
& (p%iprcparm(coarse_mat_) == repl_mat_)) )
ml_global_nmb = ( (p%iprcparm(mld_aggr_kind_) == mld_smooth_prol_).or.&
& ( (p%iprcparm(mld_aggr_kind_) == mld_biz_prol_).and.&
& (p%iprcparm(mld_coarse_mat_) == mld_repl_mat_)) )
if (ml_global_nmb) then
@ -217,9 +217,9 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
call psb_sp_scal(am3,p%dorig,info)
if(info /= 0) goto 9999
if (p%iprcparm(aggr_eig_) == max_norm_) then
if (p%iprcparm(mld_aggr_eig_) == mld_max_norm_) then
if (p%iprcparm(aggr_kind_) == biz_prol_) then
if (p%iprcparm(mld_aggr_kind_) == mld_biz_prol_) then
!
! This only works with CSR.
@ -244,15 +244,15 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
anorm = psb_spnrmi(am3,desc_a,info)
endif
omega = 4.d0/(3.d0*anorm)
p%dprcparm(aggr_damp_) = omega
p%dprcparm(mld_aggr_damp_) = omega
else if (p%iprcparm(aggr_eig_) == user_choice_) then
else if (p%iprcparm(mld_aggr_eig_) == mld_user_choice_) then
omega = p%dprcparm(aggr_damp_)
omega = p%dprcparm(mld_aggr_damp_)
else if (p%iprcparm(aggr_eig_) /= user_choice_) then
else if (p%iprcparm(mld_aggr_eig_) /= mld_user_choice_) then
write(0,*) me,'Error: invalid choice for OMEGA in blaggrmat?? ',&
& p%iprcparm(aggr_eig_)
& p%iprcparm(mld_aggr_eig_)
end if
@ -370,7 +370,7 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
call psb_numbmm(a,am1,am3)
if (debug) write(0,*) me,'Done NUMBMM 2'
if (p%iprcparm(aggr_kind_) == smooth_prol_) then
if (p%iprcparm(mld_aggr_kind_) == mld_smooth_prol_) then
call psb_transp(am1,am2,fmt='COO')
nzl = am2%infoa(psb_nnz_)
i=0
@ -398,7 +398,7 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
endif
if (debug) write(0,*) me,'starting sphalo/ rwxtd'
if (p%iprcparm(aggr_kind_) == smooth_prol_) then
if (p%iprcparm(mld_aggr_kind_) == mld_smooth_prol_) then
! am2 = ((i-wDA)Ptilde)^T
call psb_sphalo(am3,desc_a,am4,info,&
& colcnv=.false.,rowscale=.true.)
@ -418,7 +418,7 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
goto 9999
end if
else if (p%iprcparm(aggr_kind_) == biz_prol_) then
else if (p%iprcparm(mld_aggr_kind_) == mld_biz_prol_) then
call psb_rwextd(ncol,am3,info)
if(info /= 0) then
@ -454,13 +454,13 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
if (test_dump) call psb_csprt(80+me,b,head='% Smoothed aggregate AC.')
select case(p%iprcparm(aggr_kind_))
select case(p%iprcparm(mld_aggr_kind_))
case(smooth_prol_)
case(mld_smooth_prol_)
select case(p%iprcparm(coarse_mat_))
select case(p%iprcparm(mld_coarse_mat_))
case(distr_mat_)
case(mld_distr_mat_)
call psb_sp_clone(b,ac,info)
if(info /= 0) goto 9999
@ -550,7 +550,7 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
am2%m=desc_ac%matrix_data(psb_n_col_)
if (debug) write(0,*) me,'Done ac '
case(repl_mat_)
case(mld_repl_mat_)
!
!
call psb_cdrep(ntaggr,ictxt,desc_ac,info)
@ -602,11 +602,11 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
end select
case(biz_prol_)
case(mld_biz_prol_)
select case(p%iprcparm(coarse_mat_))
select case(p%iprcparm(mld_coarse_mat_))
case(distr_mat_)
case(mld_distr_mat_)
call psb_sp_clone(b,ac,info)
if(info /= 0) then
@ -632,7 +632,7 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
end if
case(repl_mat_)
case(mld_repl_mat_)
!
!

4
mlprec/mld_dasmat_bld.f90
Unescape Escape View File

@ -97,7 +97,7 @@ Subroutine mld_dasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt)
nhalo = n_col-nrow_a
If (ptype == bjac_) Then
If (ptype == mld_bjac_) Then
!
! Block Jacobi. Copy the descriptor, just in case we want to
! do the renumbering.
@ -124,7 +124,7 @@ Subroutine mld_dasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt)
end if
endif
Else If (ptype == as_) Then
Else If (ptype == mld_as_) Then
!

32
mlprec/mld_dbaseprec_aply.f90
Unescape Escape View File

@ -83,13 +83,13 @@ subroutine mld_dbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
goto 9999
end select
select case(prec%iprcparm(prec_type_))
select case(prec%iprcparm(mld_prec_type_))
case(noprec_)
case(mld_noprec_)
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
case(diag_)
case(mld_diag_)
if (size(work) >= size(x)) then
ww => work
@ -113,7 +113,7 @@ subroutine mld_dbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
end if
end if
case(bjac_)
case(mld_bjac_)
call mld_bjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
if(info.ne.0) then
@ -122,9 +122,9 @@ subroutine mld_dbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
goto 9999
end if
case(as_)
case(mld_as_)
if (prec%iprcparm(n_ovr_)==0) then
if (prec%iprcparm(mld_n_ovr_)==0) then
! shortcut: this fixes performance for RAS(0) == BJA
call mld_bjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
if(info.ne.0) then
@ -180,19 +180,19 @@ subroutine mld_dbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
tx(1:nrow_d) = x(1:nrow_d)
tx(nrow_d+1:isz) = dzero
if (prec%iprcparm(sub_restr_)==psb_halo_) then
if (prec%iprcparm(mld_sub_restr_)==psb_halo_) then
call psb_halo(tx,prec%desc_data,info,work=aux,data=psb_comm_ext_)
if(info /=0) then
info=4010
ch_err='psb_halo'
goto 9999
end if
else if (prec%iprcparm(sub_restr_) /= psb_none_) then
else if (prec%iprcparm(mld_sub_restr_) /= psb_none_) then
write(0,*) 'Problem in PRC_APLY: Unknown value for restriction ',&
&prec%iprcparm(sub_restr_)
&prec%iprcparm(mld_sub_restr_)
end if
if (prec%iprcparm(sub_ren_)>0) then
if (prec%iprcparm(mld_sub_ren_)>0) then
call psb_gelp('n',prec%perm,tx,info)
!!$ call dgelp('N',n_row,1,prec%perm,tx,isz,ww,isz,info)
if(info /=0) then
@ -209,7 +209,7 @@ subroutine mld_dbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
goto 9999
end if
if (prec%iprcparm(sub_ren_)>0) then
if (prec%iprcparm(mld_sub_ren_)>0) then
call psb_gelp('n',prec%invperm,ty,info)
!!$ call dgelp('N',n_row,1,prec%invperm,ty,isz,ww,isz,info)
if(info /=0) then
@ -219,7 +219,7 @@ subroutine mld_dbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
end if
endif
select case (prec%iprcparm(sub_prol_))
select case (prec%iprcparm(mld_sub_prol_))
case(psb_none_)
! Would work anyway, but since it's supposed to do nothing...
@ -227,7 +227,7 @@ subroutine mld_dbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
case(psb_sum_,psb_avg_)
call psb_ovrl(ty,prec%desc_data,info,&
& update=prec%iprcparm(sub_prol_),work=aux)
& update=prec%iprcparm(mld_sub_prol_),work=aux)
if(info /=0) then
info=4010
ch_err='psb_ovrl'
@ -236,7 +236,7 @@ subroutine mld_dbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
case default
write(0,*) 'Problem in PRC_APLY: Unknown value for prolongation ',&
& prec%iprcparm(sub_prol_)
& prec%iprcparm(mld_sub_prol_)
end select
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
@ -252,8 +252,8 @@ subroutine mld_dbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
endif
end if
case default
write(0,*) 'Invalid PRE%PREC ',prec%iprcparm(prec_type_),':',&
& min_prec_,noprec_,diag_,bjac_,as_
write(0,*) 'Invalid PRE%PREC ',prec%iprcparm(mld_prec_type_),':',&
& mld_min_prec_,mld_noprec_,mld_diag_,mld_bjac_,mld_as_
end select
call psb_erractionrestore(err_act)

36
mlprec/mld_dbaseprec_bld.f90
Unescape Escape View File

@ -91,14 +91,14 @@ subroutine mld_dbaseprc_bld(a,desc_a,p,info,upd)
! Should add check to ensure all procs have the same...
!
call mld_check_def(p%iprcparm(prec_type_),'base_prec',&
& diag_,is_legal_base_prec)
call mld_check_def(p%iprcparm(mld_prec_type_),'base_prec',&
& mld_diag_,is_legal_base_prec)
call psb_nullify_desc(p%desc_data)
select case(p%iprcparm(prec_type_))
case (noprec_)
select case(p%iprcparm(mld_prec_type_))
case (mld_noprec_)
! Do nothing.
call psb_cdcpy(desc_a,p%desc_data,info)
if(info /= 0) then
@ -108,7 +108,7 @@ subroutine mld_dbaseprc_bld(a,desc_a,p,info,upd)
goto 9999
end if
case (diag_)
case (mld_diag_)
call mld_diag_bld(a,desc_a,p,iupd,info)
if(debug) write(0,*)me,': out of mld_diag_bld'
@ -119,22 +119,22 @@ subroutine mld_dbaseprc_bld(a,desc_a,p,info,upd)
goto 9999
end if
case(bjac_,as_)
case(mld_bjac_,mld_as_)
call mld_check_def(p%iprcparm(n_ovr_),'overlap',&
call mld_check_def(p%iprcparm(mld_n_ovr_),'overlap',&
& 0,is_legal_n_ovr)
call mld_check_def(p%iprcparm(sub_restr_),'restriction',&
call mld_check_def(p%iprcparm(mld_sub_restr_),'restriction',&
& psb_halo_,is_legal_restrict)
call mld_check_def(p%iprcparm(sub_prol_),'prolongator',&
call mld_check_def(p%iprcparm(mld_sub_prol_),'prolongator',&
& psb_none_,is_legal_prolong)
call mld_check_def(p%iprcparm(sub_ren_),'renumbering',&
& renum_none_,is_legal_renum)
call mld_check_def(p%iprcparm(sub_solve_),'fact',&
& ilu_n_,is_legal_ml_fact)
call mld_check_def(p%iprcparm(mld_sub_ren_),'renumbering',&
& mld_renum_none_,is_legal_renum)
call mld_check_def(p%iprcparm(mld_sub_solve_),'fact',&
& mld_ilu_n_,is_legal_ml_fact)
if (p%iprcparm(sub_solve_)==sludist_) then
p%iprcparm(n_ovr_) = 0
p%iprcparm(smooth_sweeps_) = 1
if (p%iprcparm(mld_sub_solve_)==mld_sludist_) then
p%iprcparm(mld_n_ovr_) = 0
p%iprcparm(mld_smooth_sweeps_) = 1
end if
if (debug) write(0,*)me, ': Calling mld_bjac_bld'
@ -149,7 +149,7 @@ subroutine mld_dbaseprc_bld(a,desc_a,p,info,upd)
case default
info=4010
ch_err='Unknown prec_type_'
ch_err='Unknown mld_prec_type_'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
@ -157,7 +157,7 @@ subroutine mld_dbaseprc_bld(a,desc_a,p,info,upd)
p%base_a => a
p%base_desc => desc_a
p%iprcparm(prec_status_) = prec_built
p%iprcparm(mld_prec_status_) = mld_prec_built_
call psb_erractionrestore(err_act)
return

76
mlprec/mld_dbjac_aply.f90
Unescape Escape View File

@ -108,71 +108,71 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
endif
if (debug) then
write(0,*) me,' BJAC_APLY: ',prec%iprcparm(sub_solve_),prec%iprcparm(smooth_sweeps_)
write(0,*) me,' mld_bjac_APLY: ',prec%iprcparm(mld_sub_solve_),prec%iprcparm(mld_smooth_sweeps_)
end if
if (prec%iprcparm(smooth_sweeps_) == 1) then
if (prec%iprcparm(mld_smooth_sweeps_) == 1) then
select case(prec%iprcparm(sub_solve_))
case(ilu_n_,milu_n_,ilu_t_)
select case(prec%iprcparm(mld_sub_solve_))
case(mld_ilu_n_,mld_milu_n_,mld_ilu_t_)
select case(toupper(trans))
case('N')
call psb_spsm(done,prec%av(l_pr_),x,dzero,ww,desc_data,info,&
call psb_spsm(done,prec%av(mld_l_pr_),x,dzero,ww,desc_data,info,&
& trans='N',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,&
call psb_spsm(alpha,prec%av(mld_u_pr_),ww,beta,y,desc_data,info,&
& trans='N',unit='U',choice=psb_none_, work=aux)
if(info /=0) goto 9999
case('T','C')
call psb_spsm(done,prec%av(u_pr_),x,dzero,ww,desc_data,info,&
call psb_spsm(done,prec%av(mld_u_pr_),x,dzero,ww,desc_data,info,&
& 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,&
call psb_spsm(alpha,prec%av(mld_l_pr_),ww,beta,y,desc_data,info,&
& trans=trans,unit='U',choice=psb_none_,work=aux)
if(info /=0) goto 9999
end select
case(slu_)
case(mld_slu_)
ww(1:n_row) = x(1:n_row)
select case(toupper(trans))
case('N')
call mld_dslu_solve(0,n_row,1,ww,n_row,prec%iprcparm(slu_ptr_),info)
call mld_dslu_solve(0,n_row,1,ww,n_row,prec%iprcparm(mld_slu_ptr_),info)
case('T','C')
call mld_dslu_solve(1,n_row,1,ww,n_row,prec%iprcparm(slu_ptr_),info)
call mld_dslu_solve(1,n_row,1,ww,n_row,prec%iprcparm(mld_slu_ptr_),info)
end select
if(info /=0) goto 9999
call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
case(sludist_)
case(mld_sludist_)
!!$ write(0,*) 'Calling SLUDist_solve ',n_row
!!$ write(0,*) 'Calling mld_sludist_solve ',n_row
ww(1:n_row) = x(1:n_row)
select case(toupper(trans))
case('N')
call mld_dsludist_solve(0,n_row,1,ww,n_row,prec%iprcparm(slud_ptr_),info)
call mld_dsludist_solve(0,n_row,1,ww,n_row,prec%iprcparm(mld_slud_ptr_),info)
case('T','C')
call mld_dsludist_solve(1,n_row,1,ww,n_row,prec%iprcparm(slud_ptr_),info)
call mld_dsludist_solve(1,n_row,1,ww,n_row,prec%iprcparm(mld_slud_ptr_),info)
end select
if(info /=0) goto 9999
call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
case (umf_)
case (mld_umf_)
select case(toupper(trans))
case('N')
call mld_dumf_solve(0,n_row,ww,x,n_row,prec%iprcparm(umf_numptr_),info)
call mld_dumf_solve(0,n_row,ww,x,n_row,prec%iprcparm(mld_umf_numptr_),info)
case('T','C')
call mld_dumf_solve(1,n_row,ww,x,n_row,prec%iprcparm(umf_numptr_),info)
call mld_dumf_solve(1,n_row,ww,x,n_row,prec%iprcparm(mld_umf_numptr_),info)
end select
if(info /=0) goto 9999
@ -180,15 +180,15 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
case default
write(0,*) 'Unknown factorization type in bjac_aply',prec%iprcparm(sub_solve_)
write(0,*) 'Unknown factorization type in mld_bjac_aply',prec%iprcparm(mld_sub_solve_)
end select
if (debugprt) write(0,*)' Y: ',y(:)
else if (prec%iprcparm(smooth_sweeps_) > 1) then
else if (prec%iprcparm(mld_smooth_sweeps_) > 1) then
! Note: we have to add TRANS to this one !!!!!!!!!
if (size(prec%av) < ap_nd_) then
if (size(prec%av) < mld_ap_nd_) then
info = 4011
goto 9999
endif
@ -203,50 +203,50 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
tx = dzero
ty = dzero
select case(prec%iprcparm(sub_solve_))
case(ilu_n_,milu_n_,ilu_t_)
do i=1, prec%iprcparm(smooth_sweeps_)
select case(prec%iprcparm(mld_sub_solve_))
case(mld_ilu_n_,mld_milu_n_,mld_ilu_t_)
do i=1, prec%iprcparm(mld_smooth_sweeps_)
! X(k+1) = M^-1*(b-N*X(k))
ty(1:n_row) = x(1:n_row)
call psb_spmm(-done,prec%av(ap_nd_),tx,done,ty,&
call psb_spmm(-done,prec%av(mld_ap_nd_),tx,done,ty,&
& prec%desc_data,info,work=aux)
if(info /=0) goto 9999
call psb_spsm(done,prec%av(l_pr_),ty,dzero,ww,&
call psb_spsm(done,prec%av(mld_l_pr_),ty,dzero,ww,&
& prec%desc_data,info,&
& trans='N',unit='L',diag=prec%d,choice=psb_none_,work=aux)
if(info /=0) goto 9999
call psb_spsm(done,prec%av(u_pr_),ww,dzero,tx,&
call psb_spsm(done,prec%av(mld_u_pr_),ww,dzero,tx,&
& prec%desc_data,info,&
& trans='N',unit='U',choice=psb_none_,work=aux)
if(info /=0) goto 9999
end do
case(sludist_)
write(0,*) 'No sense in having SLUDist with JAC_SWEEPS >1'
case(mld_sludist_)
write(0,*) 'No sense in having SLUDist with Jmld_ac_SWEEPS >1'
info=4010
goto 9999
case(slu_)
do i=1, prec%iprcparm(smooth_sweeps_)
case(mld_slu_)
do i=1, prec%iprcparm(mld_smooth_sweeps_)
! X(k+1) = M^-1*(b-N*X(k))
ty(1:n_row) = x(1:n_row)
call psb_spmm(-done,prec%av(ap_nd_),tx,done,ty,&
call psb_spmm(-done,prec%av(mld_ap_nd_),tx,done,ty,&
& prec%desc_data,info,work=aux)
if(info /=0) goto 9999
call mld_dslu_solve(0,n_row,1,ty,n_row,prec%iprcparm(slu_ptr_),info)
call mld_dslu_solve(0,n_row,1,ty,n_row,prec%iprcparm(mld_slu_ptr_),info)
if(info /=0) goto 9999
tx(1:n_row) = ty(1:n_row)
end do
case(umf_)
do i=1, prec%iprcparm(smooth_sweeps_)
case(mld_umf_)
do i=1, prec%iprcparm(mld_smooth_sweeps_)
! X(k+1) = M^-1*(b-N*X(k))
ty(1:n_row) = x(1:n_row)
call psb_spmm(-done,prec%av(ap_nd_),tx,done,ty,&
call psb_spmm(-done,prec%av(mld_ap_nd_),tx,done,ty,&
& prec%desc_data,info,work=aux)
if(info /=0) goto 9999
call mld_dumf_solve(0,n_row,ww,ty,n_row,&
& prec%iprcparm(umf_numptr_),info)
& prec%iprcparm(mld_umf_numptr_),info)
if(info /=0) goto 9999
tx(1:n_row) = ww(1:n_row)
end do
@ -262,7 +262,7 @@ subroutine mld_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
else
info = 10
call psb_errpush(info,name,&
& i_err=(/2,prec%iprcparm(smooth_sweeps_),0,0,0/))
& i_err=(/2,prec%iprcparm(mld_smooth_sweeps_),0,0,0/))
goto 9999
endif

122
mlprec/mld_dbjac_bld.f90
Unescape Escape View File

@ -93,10 +93,10 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
endif
trans = 'N'
unitd = 'U'
if (p%iprcparm(n_ovr_) < 0) then
if (p%iprcparm(mld_n_ovr_) < 0) then
info = 11
int_err(1) = 1
int_err(2) = p%iprcparm(n_ovr_)
int_err(2) = p%iprcparm(mld_n_ovr_)
call psb_errpush(info,name,i_err=int_err)
goto 9999
endif
@ -107,9 +107,9 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
t1= psb_wtime()
if(debug) write(0,*)me,': calling mld_asmat_bld',&
& p%iprcparm(prec_type_),p%iprcparm(n_ovr_)
& p%iprcparm(mld_prec_type_),p%iprcparm(mld_n_ovr_)
if (debug) call psb_barrier(ictxt)
call mld_asmat_bld(p%iprcparm(prec_type_),p%iprcparm(n_ovr_),a,&
call mld_asmat_bld(p%iprcparm(mld_prec_type_),p%iprcparm(mld_n_ovr_),a,&
& blck,desc_a,upd,p%desc_data,info,outfmt=csrfmt)
if (debugprt) then
@ -131,7 +131,7 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
if (debug) call psb_barrier(ictxt)
select case(p%iprcparm(sub_ren_))
select case(p%iprcparm(mld_sub_ren_))
case (1:)
@ -150,23 +150,23 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
!------------------------------------------------------------------
! Split AC=M+N N off-diagonal part
! Output in COO format.
call psb_sp_clip(atmp,p%av(ap_nd_),info,&
call psb_sp_clip(atmp,p%av(mld_ap_nd_),info,&
& jmin=atmp%m+1,rscale=.false.,cscale=.false.)
call psb_spcnv(p%av(ap_nd_),info,afmt='csr',dupl=psb_dupl_add_)
call psb_spcnv(p%av(mld_ap_nd_),info,afmt='csr',dupl=psb_dupl_add_)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_spcnv csr 1')
goto 9999
end if
k = psb_sp_get_nnzeros(p%av(ap_nd_))
k = psb_sp_get_nnzeros(p%av(mld_ap_nd_))
call psb_sum(ictxt,k)
if (k == 0) then
! If the off diagonal part is emtpy, there's no point
! in doing multiple Jacobi sweeps. This is certain
! to happen when running on a single processor.
p%iprcparm(smooth_sweeps_) = 1
p%iprcparm(mld_smooth_sweeps_) = 1
end if
@ -180,9 +180,9 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
if (debug) write(0,*) me,' Factoring rows ',&
&atmp%m,a%m,blck%m,atmp%ia2(atmp%m+1)-1
select case(p%iprcparm(sub_solve_))
select case(p%iprcparm(mld_sub_solve_))
case(ilu_n_,milu_n_,ilu_t_)
case(mld_ilu_n_,mld_milu_n_,mld_ilu_t_)
call psb_spcnv(atmp,info,afmt='csr',dupl=psb_dupl_add_)
if (info /= 0) then
@ -202,18 +202,18 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
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
call psb_csprt(80+me,p%av(mld_l_pr_),head='% Local L factor')
write(80+me,*) '% Diagonal: ',p%av(mld_l_pr_)%m
do i=1,p%av(mld_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')
call psb_csprt(80+me,p%av(mld_u_pr_),head='% Local U factor')
close(80+me)
endif
case(slu_)
case(mld_slu_)
call psb_spcnv(atmp,info,afmt='csr',dupl=psb_dupl_add_)
if (info /= 0) then
@ -223,11 +223,11 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
call mld_slu_bld(atmp,p%desc_data,p,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='slu_bld')
call psb_errpush(4010,name,a_err='mld_slu_bld')
goto 9999
end if
case(umf_)
case(mld_umf_)
call psb_spcnv(atmp,info,afmt='csc',dupl=psb_dupl_add_)
if (info /= 0) then
@ -236,20 +236,20 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
end if
call mld_umf_bld(atmp,p%desc_data,p,info)
if(debug) write(0,*)me,': Done umf_bld ',info
if(debug) write(0,*)me,': Done mld_umf_bld ',info
if (info /= 0) then
call psb_errpush(4010,name,a_err='umf_bld')
call psb_errpush(4010,name,a_err='mld_umf_bld')
goto 9999
end if
case(f_none_)
case(mld_f_none_)
info=4010
call psb_errpush(info,name,a_err='Inconsistent prec f_none_')
call psb_errpush(info,name,a_err='Inconsistent prec mld_f_none_')
goto 9999
case default
info=4010
call psb_errpush(info,name,a_err='Unknown sub_solve_')
call psb_errpush(info,name,a_err='Unknown mld_sub_solve_')
goto 9999
end select
@ -266,37 +266,37 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
case(0) ! No renumbering
select case(p%iprcparm(sub_solve_))
select case(p%iprcparm(mld_sub_solve_))
case(ilu_n_,milu_n_,ilu_t_)
case(mld_ilu_n_,mld_milu_n_,mld_ilu_t_)
if (p%iprcparm(smooth_sweeps_) > 1) then
if (p%iprcparm(mld_smooth_sweeps_) > 1) then
n_row = psb_cd_get_local_rows(p%desc_data)
n_col = psb_cd_get_local_cols(p%desc_data)
nrow_a = a%m
! The following is known to work
! given that the output from CLIP is in COO.
call psb_sp_clip(a,p%av(ap_nd_),info,&
call psb_sp_clip(a,p%av(mld_ap_nd_),info,&
& jmin=nrow_a+1,rscale=.false.,cscale=.false.)
call psb_sp_clip(blck,atmp,info,&
& jmin=nrow_a+1,rscale=.false.,cscale=.false.)
call psb_rwextd(n_row,p%av(ap_nd_),info,b=atmp)
call psb_rwextd(n_row,p%av(mld_ap_nd_),info,b=atmp)
call psb_spcnv(p%av(ap_nd_),info,afmt='csr',dupl=psb_dupl_add_)
call psb_spcnv(p%av(mld_ap_nd_),info,afmt='csr',dupl=psb_dupl_add_)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_spcnv csr 4')
goto 9999
end if
k = psb_sp_get_nnzeros(p%av(ap_nd_))
k = psb_sp_get_nnzeros(p%av(mld_ap_nd_))
call psb_sum(ictxt,k)
if (k == 0) then
! If the off diagonal part is emtpy, there's no point
! in doing multiple Jacobi sweeps. This is certain
! to happen when running on a single processor.
p%iprcparm(smooth_sweeps_) = 1
p%iprcparm(mld_smooth_sweeps_) = 1
end if
call psb_sp_free(atmp,info)
end if
@ -313,18 +313,18 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
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
call psb_csprt(80+me,p%av(mld_l_pr_),head='% Local L factor')
write(80+me,*) '% Diagonal: ',p%av(mld_l_pr_)%m
do i=1,p%av(mld_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')
call psb_csprt(80+me,p%av(mld_u_pr_),head='% Local U factor')
close(80+me)
endif
case(slu_)
case(mld_slu_)
call psb_spcnv(a,atmp,info,afmt='coo')
if (info /= 0) then
@ -336,34 +336,34 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
n_col = psb_cd_get_local_cols(p%desc_data)
call psb_rwextd(n_row,atmp,info,b=blck)
if (p%iprcparm(smooth_sweeps_) > 1) then
if (p%iprcparm(mld_smooth_sweeps_) > 1) then
!------------------------------------------------------------------
! Split AC=M+N N off-diagonal part
! Output in COO format.
call psb_sp_clip(atmp,p%av(ap_nd_),info,&
call psb_sp_clip(atmp,p%av(mld_ap_nd_),info,&
& jmin=atmp%m+1,rscale=.false.,cscale=.false.)
call psb_spcnv(p%av(ap_nd_),info,afmt='csr',dupl=psb_dupl_add_)
call psb_spcnv(p%av(mld_ap_nd_),info,afmt='csr',dupl=psb_dupl_add_)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_spcnv csr 6')
goto 9999
end if
k = psb_sp_get_nnzeros(p%av(ap_nd_))
k = psb_sp_get_nnzeros(p%av(mld_ap_nd_))
call psb_sum(ictxt,k)
if (k == 0) then
! If the off diagonal part is emtpy, there's no point
! in doing multiple Jacobi sweeps. This is certain
! to happen when running on a single processor.
p%iprcparm(smooth_sweeps_) = 1
p%iprcparm(mld_smooth_sweeps_) = 1
end if
endif
if (info == 0) call psb_spcnv(atmp,info,afmt='csr',dupl=psb_dupl_add_)
if (info == 0) call mld_slu_bld(atmp,p%desc_data,p,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='slu_bld')
call psb_errpush(4010,name,a_err='mld_slu_bld')
goto 9999
end if
@ -374,7 +374,7 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
end if
case(sludist_)
case(mld_sludist_)
call psb_spcnv(a,atmp,info,afmt='coo')
if (info /= 0) then
@ -386,34 +386,34 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
n_col = psb_cd_get_local_cols(p%desc_data)
call psb_rwextd(n_row,atmp,info,b=blck)
if (p%iprcparm(smooth_sweeps_) > 1) then
if (p%iprcparm(mld_smooth_sweeps_) > 1) then
!------------------------------------------------------------------
! Split AC=M+N N off-diagonal part
! Output in COO format.
call psb_sp_clip(atmp,p%av(ap_nd_),info,&
call psb_sp_clip(atmp,p%av(mld_ap_nd_),info,&
& jmin=atmp%m+1,rscale=.false.,cscale=.false.)
call psb_spcnv(p%av(ap_nd_),info,afmt='csr',dupl=psb_dupl_add_)
call psb_spcnv(p%av(mld_ap_nd_),info,afmt='csr',dupl=psb_dupl_add_)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_spcnv csr 7')
goto 9999
end if
k = psb_sp_get_nnzeros(p%av(ap_nd_))
k = psb_sp_get_nnzeros(p%av(mld_ap_nd_))
call psb_sum(ictxt,k)
if (k == 0) then
! If the off diagonal part is emtpy, there's no point
! in doing multiple Jacobi sweeps. This is certain
! to happen when running on a single processor.
p%iprcparm(smooth_sweeps_) = 1
p%iprcparm(mld_smooth_sweeps_) = 1
end if
endif
if (info == 0) call psb_spcnv(atmp,info,afmt='csr',dupl=psb_dupl_add_)
if (info == 0) call mld_sludist_bld(atmp,p%desc_data,p,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='slu_bld')
call psb_errpush(4010,name,a_err='mld_slu_bld')
goto 9999
end if
@ -423,7 +423,7 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
goto 9999
end if
case(umf_)
case(mld_umf_)
call psb_spcnv(a,atmp,info,afmt='coo')
@ -436,28 +436,28 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
n_col = psb_cd_get_local_cols(p%desc_data)
call psb_rwextd(n_row,atmp,info,b=blck)
if (p%iprcparm(smooth_sweeps_) > 1) then
if (p%iprcparm(mld_smooth_sweeps_) > 1) then
!------------------------------------------------------------------
! Split AC=M+N N off-diagonal part
! Output in COO format.
!!$ write(0,*) 'bjac_bld:' size(p%av),ap_nd_
call psb_sp_clip(atmp,p%av(ap_nd_),info,&
!!$ write(0,*) 'mld_bjac_bld:' size(p%av),mld_ap_nd_
call psb_sp_clip(atmp,p%av(mld_ap_nd_),info,&
& jmin=atmp%m+1,rscale=.false.,cscale=.false.)
call psb_spcnv(p%av(ap_nd_),info,afmt='csr',dupl=psb_dupl_add_)
call psb_spcnv(p%av(mld_ap_nd_),info,afmt='csr',dupl=psb_dupl_add_)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_spcnv csr 8')
goto 9999
end if
k = psb_sp_get_nnzeros(p%av(ap_nd_))
k = psb_sp_get_nnzeros(p%av(mld_ap_nd_))
call psb_sum(ictxt,k)
if (k == 0) then
! If the off diagonal part is emtpy, there's no point
! in doing multiple Jacobi sweeps. This is certain
! to happen when running on a single processor.
p%iprcparm(smooth_sweeps_) = 1
p%iprcparm(mld_smooth_sweeps_) = 1
end if
endif
@ -468,9 +468,9 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
end if
call mld_umf_bld(atmp,p%desc_data,p,info)
if(debug) write(0,*)me,': Done umf_bld ',info
if(debug) write(0,*)me,': Done mld_umf_bld ',info
if (info /= 0) then
call psb_errpush(4010,name,a_err='umf_bld')
call psb_errpush(4010,name,a_err='mld_umf_bld')
goto 9999
end if
@ -481,14 +481,14 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info)
end if
case(f_none_)
case(mld_f_none_)
info=4010
call psb_errpush(info,name,a_err='Inconsistent prec f_none_')
call psb_errpush(info,name,a_err='Inconsistent prec mld_f_none_')
goto 9999
case default
info=4010
call psb_errpush(info,name,a_err='Unknown sub_solve_')
call psb_errpush(info,name,a_err='Unknown mld_sub_solve_')
goto 9999
end select

34
mlprec/mld_dilu_bld.f90
Unescape Escape View File

@ -87,7 +87,7 @@ subroutine mld_dilu_bld(a,desc_a,p,upd,info,blck)
unitd = 'U'
if (allocated(p%av)) then
if (size(p%av) < bp_ilu_avsz) then
if (size(p%av) < mld_bp_ilu_avsz_) then
do i=1,size(p%av)
call psb_sp_free(p%av(i),info)
if (info /= 0) then
@ -100,7 +100,7 @@ subroutine mld_dilu_bld(a,desc_a,p,upd,info,blck)
endif
end if
if (.not.allocated(p%av)) then
allocate(p%av(max_avsz),stat=info)
allocate(p%av(mld_max_avsz_),stat=info)
if (info /= 0) then
call psb_errpush(4000,name)
goto 9999
@ -117,12 +117,12 @@ subroutine mld_dilu_bld(a,desc_a,p,upd,info,blck)
if (debug) call psb_barrier(ictxt)
n_row = p%desc_data%matrix_data(psb_n_row_)
p%av(l_pr_)%m = n_row
p%av(l_pr_)%k = n_row
p%av(u_pr_)%m = n_row
p%av(u_pr_)%k = n_row
call psb_sp_all(n_row,n_row,p%av(l_pr_),nztota,info)
if (info == 0) call psb_sp_all(n_row,n_row,p%av(u_pr_),nztota,info)
p%av(mld_l_pr_)%m = n_row
p%av(mld_l_pr_)%k = n_row
p%av(mld_u_pr_)%m = n_row
p%av(mld_u_pr_)%k = n_row
call psb_sp_all(n_row,n_row,p%av(mld_l_pr_),nztota,info)
if (info == 0) call psb_sp_all(n_row,n_row,p%av(mld_u_pr_),nztota,info)
if(info/=0) then
info=4010
ch_err='psb_sp_all'
@ -149,7 +149,7 @@ subroutine mld_dilu_bld(a,desc_a,p,upd,info,blck)
! Ok, factor the matrix.
!
t5 = psb_wtime()
call mld_ilu_fct(p%iprcparm(sub_solve_),a,p%av(l_pr_),p%av(u_pr_),&
call mld_ilu_fct(p%iprcparm(mld_sub_solve_),a,p%av(mld_l_pr_),p%av(mld_u_pr_),&
& p%d,info,blck=blck)
if(info/=0) then
info=4010
@ -165,12 +165,12 @@ subroutine mld_dilu_bld(a,desc_a,p,upd,info,blck)
!
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
call psb_csprt(80+me,p%av(mld_l_pr_),head='% Local L factor')
write(80+me,*) '% Diagonal: ',p%av(mld_l_pr_)%m
do i=1,p%av(mld_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')
call psb_csprt(80+me,p%av(mld_u_pr_),head='% Local U factor')
close(80+me)
endif
@ -184,12 +184,12 @@ subroutine mld_dilu_bld(a,desc_a,p,upd,info,blck)
! 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)
if (psb_sp_getifld(psb_upd_,p%av(mld_u_pr_),info) /= psb_upd_perm_) then
call psb_sp_trim(p%av(mld_u_pr_),info)
endif
if (psb_sp_getifld(psb_upd_,p%av(l_pr_),info) /= psb_upd_perm_) then
call psb_sp_trim(p%av(l_pr_),info)
if (psb_sp_getifld(psb_upd_,p%av(mld_l_pr_),info) /= psb_upd_perm_) then
call psb_sp_trim(p%av(mld_l_pr_),info)
endif

4
mlprec/mld_dilu_fct.f90
Unescape Escape View File

@ -291,7 +291,7 @@ contains
! which means that this entry does not match; thus
! we take it out of diagonal for MILU.
!
if (ialg == milu_n_) then
if (ialg == mld_milu_n_) then
dia = dia - temp*uaspk(jj)
end if
enddo updateloop
@ -435,7 +435,7 @@ contains
! which means that this entry does not match; thus
! we take it out of diagonal for MILU.
!
if (ialg == milu_n_) then
if (ialg == mld_milu_n_) then
dia = dia - temp*uaspk(jj)
end if
enddo updateloopb

116
mlprec/mld_dmlprec_aply.f90
Unescape Escape View File

@ -61,11 +61,11 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
! 1. Number of levels = NLEV = size(baseprecv(:))
! 2. baseprecv(ilev)%av(:) sparse matrices needed for the current level.
! Includes:
! 2.1.: baseprecv(ilev)%av(l_pr_) L factor of ILU preconditioners
! 2.2.: baseprecv(ilev)%av(u_pr_) U factor of ILU preconditioners
! 2.3.: baseprecv(ilev)%av(ap_nd_) Off-diagonal part of A for Jacobi sweeps
! 2.4.: baseprecv(ilev)%av(ac_) Aggregated matrix of level ILEV
! 2.5.: baseprecv(ilev)%av(sm_pr_t_) Smoother prolongator transpose; maps vectors
! 2.1.: baseprecv(ilev)%av(mld_l_pr_) L factor of ILU preconditioners
! 2.2.: baseprecv(ilev)%av(mld_u_pr_) U factor of ILU preconditioners
! 2.3.: baseprecv(ilev)%av(mld_ap_nd_) Off-diagonal part of A for Jacobi sweeps
! 2.4.: baseprecv(ilev)%av(mld_ac_) Aggregated matrix of level ILEV
! 2.5.: baseprecv(ilev)%av(mld_sm_pr_t_) Smoother prolongator transpose; maps vectors
! (ilev-1) ---> (ilev)
! 2.6.: baseprecv(ilev)%av(sm_pr_) Smoother prolongator; maps vectors
! (ilev) ---> (ilev-1)
@ -74,7 +74,7 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
! 3. baseprecv(ilev)%desc_data comm descriptor for level ILEV
! 4. baseprecv(ilev)%base_a Pointer (really a pointer!) to the base matrix
! baseprecv(ilev)%base_desc of the current level, i.e.: if ILEV=1 then A
! else the aggregated matrix av(ac_); so we have
! else the aggregated matrix av(mld_ac_); so we have
! a unified treatment of residuals. Need this to
! avoid passing explicitly matrix A to the
! outer prec. routine
@ -130,15 +130,15 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
end if
select case(baseprecv(2)%iprcparm(ml_type_))
select case(baseprecv(2)%iprcparm(mld_ml_type_))
case(no_ml_)
case(mld_no_ml_)
! Should not really get here.
call psb_errpush(4010,name,a_err='no_ml_ in mlprc_aply?')
call psb_errpush(4010,name,a_err='mld_no_ml_ in mlprc_aply?')
goto 9999
case(add_ml_)
case(mld_add_ml_)
!
@ -190,9 +190,9 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
mlprec_wrk(ilev)%tx(n_row+1:max(n_row,n_col)) = dzero
mlprec_wrk(ilev)%ty(:) = dzero
ismth = baseprecv(ilev)%iprcparm(aggr_kind_)
icm = baseprecv(ilev)%iprcparm(coarse_mat_)
if (ismth /= no_smooth_) then
ismth = baseprecv(ilev)%iprcparm(mld_aggr_kind_)
icm = baseprecv(ilev)%iprcparm(mld_coarse_mat_)
if (ismth /= mld_no_smooth_) then
!
! Smoothed aggregation
!
@ -200,7 +200,7 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
& info,work=work)
if(info /=0) goto 9999
call psb_csmm(done,baseprecv(ilev)%av(sm_pr_t_),mlprec_wrk(ilev-1)%x2l,&
call psb_csmm(done,baseprecv(ilev)%av(mld_sm_pr_t_),mlprec_wrk(ilev-1)%x2l,&
& dzero,mlprec_wrk(ilev)%x2l,info)
if(info /=0) goto 9999
@ -216,10 +216,10 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
end if
if (icm ==repl_mat_) Then
if (icm ==mld_repl_mat_) Then
call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nr2l))
else if (icm/= distr_mat_) Then
write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',icm
else if (icm/= mld_distr_mat_) Then
write(0,*) 'Unknown value for baseprecv(2)%iprcparm(mld_coarse_mat_) ',icm
endif
call mld_baseprec_aply(done,baseprecv(ilev),&
@ -234,10 +234,10 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
n_col = psb_cd_get_local_cols(baseprecv(ilev-1)%desc_data)
nc2l = psb_cd_get_local_cols(baseprecv(ilev)%desc_data)
nr2l = psb_cd_get_local_rows(baseprecv(ilev)%desc_data)
ismth = baseprecv(ilev)%iprcparm(aggr_kind_)
icm = baseprecv(ilev)%iprcparm(coarse_mat_)
ismth = baseprecv(ilev)%iprcparm(mld_aggr_kind_)
icm = baseprecv(ilev)%iprcparm(mld_coarse_mat_)
if (ismth /= no_smooth_) then
if (ismth /= mld_no_smooth_) then
call psb_csmm(done,baseprecv(ilev)%av(sm_pr_),mlprec_wrk(ilev)%y2l,&
& done,mlprec_wrk(ilev-1)%y2l,info)
@ -257,16 +257,16 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
if(info /=0) goto 9999
case(mult_ml)
case(mld_mult_ml_)
!
! Multiplicative multilevel
! Pre/post smoothing versions.
!
select case(baseprecv(2)%iprcparm(smooth_pos_))
select case(baseprecv(2)%iprcparm(mld_smooth_pos_))
case(post_smooth_)
case(mld_post_smooth_)
!
! Post smoothing.
@ -309,8 +309,8 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
n_col = psb_cd_get_local_cols(baseprecv(ilev-1)%desc_data)
nc2l = psb_cd_get_local_cols(baseprecv(ilev)%desc_data)
nr2l = psb_cd_get_local_rows(baseprecv(ilev)%desc_data)
ismth = baseprecv(ilev)%iprcparm(aggr_kind_)
icm = baseprecv(ilev)%iprcparm(coarse_mat_)
ismth = baseprecv(ilev)%iprcparm(mld_aggr_kind_)
icm = baseprecv(ilev)%iprcparm(mld_coarse_mat_)
if (debug) write(0,*) me, 'mlpr_aply starting up sweep ',&
& ilev,allocated(baseprecv(ilev)%iprcparm),n_row,n_col,&
@ -329,7 +329,7 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
mlprec_wrk(ilev)%x2l(:) = dzero
mlprec_wrk(ilev)%y2l(:) = dzero
mlprec_wrk(ilev)%tx(:) = dzero
if (ismth /= no_smooth_) then
if (ismth /= mld_no_smooth_) then
!
! Smoothed aggregation
!
@ -340,7 +340,7 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
if(info /=0) goto 9999
if (debug) write(0,*) me, 'mlpr_aply csmm in up sweep ', ilev
call psb_csmm(done,baseprecv(ilev)%av(sm_pr_t_),mlprec_wrk(ilev-1)%x2l, &
call psb_csmm(done,baseprecv(ilev)%av(mld_sm_pr_t_),mlprec_wrk(ilev-1)%x2l, &
& dzero,mlprec_wrk(ilev)%x2l,info)
if(info /=0) goto 9999
@ -356,18 +356,18 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
end if
if (debug) write(0,*) me, 'mlpr_aply possible sum in up sweep ', &
& ilev,icm,associated(baseprecv(ilev)%base_desc),repl_mat_
& ilev,icm,associated(baseprecv(ilev)%base_desc),mld_repl_mat_
if (debug) write(0,*) me, 'mlpr_aply geaxpby in up sweep X', &
& ilev,associated(baseprecv(ilev)%base_desc),&
& baseprecv(ilev)%base_desc%matrix_data(psb_n_row_),&
& baseprecv(ilev)%base_desc%matrix_data(psb_n_col_),&
& size(mlprec_wrk(ilev)%tx),size(mlprec_wrk(ilev)%x2l)
if (icm == repl_mat_) Then
if (icm == mld_repl_mat_) Then
if (debug) write(0,*) 'Entering psb_sum ',nr2l
call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nr2l))
else if (icm /= distr_mat_) Then
write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ', icm
else if (icm /= mld_distr_mat_) Then
write(0,*) 'Unknown value for baseprecv(2)%iprcparm(mld_coarse_mat_) ', icm
endif
call psb_geaxpby(done,mlprec_wrk(ilev)%x2l,dzero,mlprec_wrk(ilev)%tx,&
@ -389,11 +389,11 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
do ilev=nlev-1, 1, -1
if (debug) write(0,*) me, 'mlpr_aply starting down sweep',ilev
ismth = baseprecv(ilev+1)%iprcparm(aggr_kind_)
ismth = baseprecv(ilev+1)%iprcparm(mld_aggr_kind_)
n_row = psb_cd_get_local_rows(baseprecv(ilev)%base_desc)
if (ismth /= no_smooth_) then
if (ismth == smooth_prol_) &
if (ismth /= mld_no_smooth_) then
if (ismth == mld_smooth_prol_) &
& call psb_halo(mlprec_wrk(ilev+1)%y2l,baseprecv(ilev+1)%desc_data,&
& info,work=work)
call psb_csmm(done,baseprecv(ilev+1)%av(sm_pr_),mlprec_wrk(ilev+1)%y2l,&
@ -426,7 +426,7 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
if(info /=0) goto 9999
case(pre_smooth_)
case(mld_pre_smooth_)
!
! Pre smoothing.
@ -479,8 +479,8 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
n_col = psb_cd_get_local_cols(baseprecv(ilev-1)%desc_data)
nc2l = psb_cd_get_local_cols(baseprecv(ilev)%desc_data)
nr2l = psb_cd_get_local_rows(baseprecv(ilev)%desc_data)
ismth = baseprecv(ilev)%iprcparm(aggr_kind_)
icm = baseprecv(ilev)%iprcparm(coarse_mat_)
ismth = baseprecv(ilev)%iprcparm(mld_aggr_kind_)
icm = baseprecv(ilev)%iprcparm(mld_coarse_mat_)
allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%y2l(nc2l),&
& mlprec_wrk(ilev)%x2l(nc2l), stat=info)
@ -496,7 +496,7 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
mlprec_wrk(ilev)%tx(:) = dzero
if (ismth /= no_smooth_) then
if (ismth /= mld_no_smooth_) then
!
!Smoothed Aggregation
!
@ -504,7 +504,7 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
& info,work=work)
if(info /=0) goto 9999
call psb_csmm(done,baseprecv(ilev)%av(sm_pr_t_),mlprec_wrk(ilev-1)%tx,dzero,&
call psb_csmm(done,baseprecv(ilev)%av(mld_sm_pr_t_),mlprec_wrk(ilev-1)%tx,dzero,&
& mlprec_wrk(ilev)%x2l,info)
if(info /=0) goto 9999
@ -520,10 +520,10 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
end do
end if
if (icm ==repl_mat_) then
if (icm ==mld_repl_mat_) then
call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nr2l))
else if (icm /= distr_mat_) then
write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ', icm
else if (icm /= mld_distr_mat_) then
write(0,*) 'Unknown value for baseprecv(2)%iprcparm(mld_coarse_mat_) ', icm
endif
@ -543,12 +543,12 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
do ilev = nlev-1, 1, -1
ismth = baseprecv(ilev+1)%iprcparm(aggr_kind_)
ismth = baseprecv(ilev+1)%iprcparm(mld_aggr_kind_)
n_row = psb_cd_get_local_rows(baseprecv(ilev)%base_desc)
if (ismth /= no_smooth_) then
if (ismth /= mld_no_smooth_) then
if (ismth == smooth_prol_) &
if (ismth == mld_smooth_prol_) &
& call psb_halo(mlprec_wrk(ilev+1)%y2l,&
& baseprecv(ilev+1)%desc_data,info,work=work)
call psb_csmm(done,baseprecv(ilev+1)%av(sm_pr_),mlprec_wrk(ilev+1)%y2l,&
@ -574,7 +574,7 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
case(twoside_smooth_)
case(mld_twoside_smooth_)
!
! Symmetrized smoothing.
@ -635,8 +635,8 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
n_col = psb_cd_get_local_cols(baseprecv(ilev-1)%desc_data)
nc2l = psb_cd_get_local_cols(baseprecv(ilev)%desc_data)
nr2l = psb_cd_get_local_rows(baseprecv(ilev)%desc_data)
ismth = baseprecv(ilev)%iprcparm(aggr_kind_)
icm = baseprecv(ilev)%iprcparm(coarse_mat_)
ismth = baseprecv(ilev)%iprcparm(mld_aggr_kind_)
icm = baseprecv(ilev)%iprcparm(mld_coarse_mat_)
allocate(mlprec_wrk(ilev)%ty(nc2l),mlprec_wrk(ilev)%y2l(nc2l),&
& mlprec_wrk(ilev)%x2l(nc2l), stat=info)
@ -653,7 +653,7 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
mlprec_wrk(ilev)%ty(:) = dzero
if (ismth /= no_smooth_) then
if (ismth /= mld_no_smooth_) then
!
!Smoothed Aggregation
!
@ -661,7 +661,7 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
& info,work=work)
if(info /=0) goto 9999
call psb_csmm(done,baseprecv(ilev)%av(sm_pr_t_),mlprec_wrk(ilev-1)%ty,dzero,&
call psb_csmm(done,baseprecv(ilev)%av(mld_sm_pr_t_),mlprec_wrk(ilev-1)%ty,dzero,&
& mlprec_wrk(ilev)%x2l,info)
if(info /=0) goto 9999
@ -677,10 +677,10 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
end do
end if
if (icm == repl_mat_) then
if (icm == mld_repl_mat_) then
call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nr2l))
else if (icm /= distr_mat_) then
write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ', icm
else if (icm /= mld_distr_mat_) then
write(0,*) 'Unknown value for baseprecv(2)%iprcparm(mld_coarse_mat_) ', icm
endif
call psb_geaxpby(done,mlprec_wrk(ilev)%x2l,dzero,mlprec_wrk(ilev)%tx,&
@ -704,11 +704,11 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
do ilev=nlev-1, 1, -1
ismth = baseprecv(ilev+1)%iprcparm(aggr_kind_)
ismth = baseprecv(ilev+1)%iprcparm(mld_aggr_kind_)
n_row = psb_cd_get_local_rows(baseprecv(ilev)%base_desc)
if (ismth /= no_smooth_) then
if (ismth == smooth_prol_) &
if (ismth /= mld_no_smooth_) then
if (ismth == mld_smooth_prol_) &
& call psb_halo(mlprec_wrk(ilev+1)%y2l,baseprecv(ilev+1)%desc_data,&
& info,work=work)
call psb_csmm(done,baseprecv(ilev+1)%av(sm_pr_),mlprec_wrk(ilev+1)%y2l,&
@ -743,14 +743,14 @@ subroutine mld_dmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
case default
call psb_errpush(4013,name,a_err='wrong smooth_pos',&
& i_Err=(/baseprecv(2)%iprcparm(smooth_pos_),0,0,0,0/))
& i_Err=(/baseprecv(2)%iprcparm(mld_smooth_pos_),0,0,0,0/))
goto 9999
end select
case default
call psb_errpush(4013,name,a_err='wrong mltype',&
& i_Err=(/baseprecv(2)%iprcparm(ml_type_),0,0,0,0/))
& i_Err=(/baseprecv(2)%iprcparm(mld_ml_type_),0,0,0,0/))
goto 9999
end select

42
mlprec/mld_dmlprec_bld.f90
Unescape Escape View File

@ -67,34 +67,34 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info)
call psb_errpush(info,name)
goto 9999
endif
call mld_check_def(p%iprcparm(ml_type_),'Multilevel type',&
& mult_ml,is_legal_ml_type)
call mld_check_def(p%iprcparm(aggr_alg_),'aggregation',&
& dec_aggr_,is_legal_ml_aggr_kind)
call mld_check_def(p%iprcparm(aggr_kind_),'Smoother kind',&
& smooth_prol_,is_legal_ml_smth_kind)
call mld_check_def(p%iprcparm(coarse_mat_),'Coarse matrix',&
& distr_mat_,is_legal_ml_coarse_mat)
call mld_check_def(p%iprcparm(smooth_pos_),'smooth_pos',&
& pre_smooth_,is_legal_ml_smooth_pos)
call mld_check_def(p%iprcparm(mld_ml_type_),'Multilevel type',&
& mld_mult_ml_,is_legal_ml_type)
call mld_check_def(p%iprcparm(mld_aggr_alg_),'aggregation',&
& mld_dec_aggr_,is_legal_ml_aggr_kind)
call mld_check_def(p%iprcparm(mld_aggr_kind_),'Smoother kind',&
& mld_smooth_prol_,is_legal_ml_smth_kind)
call mld_check_def(p%iprcparm(mld_coarse_mat_),'Coarse matrix',&
& mld_distr_mat_,is_legal_ml_coarse_mat)
call mld_check_def(p%iprcparm(mld_smooth_pos_),'smooth_pos',&
& mld_pre_smooth_,is_legal_ml_smooth_pos)
!!$ nullify(p%desc_data)
select case(p%iprcparm(sub_solve_))
case(ilu_n_)
call mld_check_def(p%iprcparm(sub_fill_in_),'Level',0,is_legal_ml_lev)
case(ilu_t_)
call mld_check_def(p%dprcparm(fact_eps_),'Eps',dzero,is_legal_ml_eps)
select case(p%iprcparm(mld_sub_solve_))
case(mld_ilu_n_)
call mld_check_def(p%iprcparm(mld_sub_fill_in_),'Level',0,is_legal_ml_lev)
case(mld_ilu_t_)
call mld_check_def(p%dprcparm(mld_fact_eps_),'Eps',dzero,is_legal_ml_eps)
end select
call mld_check_def(p%dprcparm(aggr_damp_),'omega',dzero,is_legal_omega)
call mld_check_def(p%iprcparm(smooth_sweeps_),'Jacobi sweeps',&
call mld_check_def(p%dprcparm(mld_aggr_damp_),'omega',dzero,is_legal_omega)
call mld_check_def(p%iprcparm(mld_smooth_sweeps_),'Jacobi sweeps',&
& 1,is_legal_jac_sweeps)
! Currently this is ignored by gen_aggrmap, but it could be
! changed in the future. Need to package nlaggr & mlia in a
! private data structure?
call mld_aggrmap_bld(p%iprcparm(aggr_alg_),a,desc_a,p%nlaggr,p%mlia,info)
call mld_aggrmap_bld(p%iprcparm(mld_aggr_alg_),a,desc_a,p%nlaggr,p%mlia,info)
if(info /= 0) then
info=4010
ch_err='psb_gen_aggrmap'
@ -130,11 +130,11 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info)
! We have used a separate ac because:
! 1. We want to reuse the same routines mld_ilu_bld etc.
! 2. We do NOT want to pass an argument twice to them
! p%av(ac_) and p, as this would violate the Fortran standard
! p%av(mld_ac_) and p, as this would violate the Fortran standard
! Hence a separate AC and a TRANSFER function at the end.
!
call psb_sp_transfer(ac,p%av(ac_),info)
p%base_a => p%av(ac_)
call psb_sp_transfer(ac,p%av(mld_ac_),info)
p%base_a => p%av(mld_ac_)
call psb_cdtransfer(desc_ac,p%desc_ac,info)
if (info /= 0) then

2
mlprec/mld_dprecbld.f90
Unescape Escape View File

@ -163,7 +163,7 @@ contains
if (allocated(p%av)) then
! Have not decided what to do yet
end if
allocate(p%av(max_avsz),stat=info)
allocate(p%av(mld_max_avsz_),stat=info)
!!$ if (info /= 0) return
do k=1,size(p%av)
call psb_nullify_sp(p%av(k))

164
mlprec/mld_dprecinit.f90
Unescape Escape View File

@ -61,67 +61,67 @@ subroutine mld_dprecinit(p,ptype,info,nlev)
nlev_ = 1
ilev_ = 1
allocate(p%baseprecv(nlev_),stat=info)
if (info == 0) call psb_realloc(ifpsz,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(dfpsz,p%baseprecv(ilev_)%dprcparm,info)
if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(mld_dfpsz_,p%baseprecv(ilev_)%dprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(prec_type_) = noprec_
p%baseprecv(ilev_)%iprcparm(sub_solve_) = f_none_
p%baseprecv(ilev_)%iprcparm(sub_restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(sub_prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(sub_ren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(smooth_sweeps_) = 1
p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_noprec_
p%baseprecv(ilev_)%iprcparm(mld_sub_solve_) = mld_f_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_ren_) = 0
p%baseprecv(ilev_)%iprcparm(mld_n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(mld_smooth_sweeps_) = 1
case ('DIAG')
nlev_ = 1
ilev_ = 1
allocate(p%baseprecv(nlev_),stat=info)
if (info == 0) call psb_realloc(ifpsz,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(dfpsz,p%baseprecv(ilev_)%dprcparm,info)
if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(mld_dfpsz_,p%baseprecv(ilev_)%dprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(prec_type_) = diag_
p%baseprecv(ilev_)%iprcparm(sub_solve_) = f_none_
p%baseprecv(ilev_)%iprcparm(sub_restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(sub_prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(sub_ren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(smooth_sweeps_) = 1
p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_diag_
p%baseprecv(ilev_)%iprcparm(mld_sub_solve_) = mld_f_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_ren_) = 0
p%baseprecv(ilev_)%iprcparm(mld_n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(mld_smooth_sweeps_) = 1
case ('BJAC')
nlev_ = 1
ilev_ = 1
allocate(p%baseprecv(nlev_),stat=info)
if (info == 0) call psb_realloc(ifpsz,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(dfpsz,p%baseprecv(ilev_)%dprcparm,info)
if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(mld_dfpsz_,p%baseprecv(ilev_)%dprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(prec_type_) = bjac_
p%baseprecv(ilev_)%iprcparm(sub_solve_) = ilu_n_
p%baseprecv(ilev_)%iprcparm(sub_restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(sub_prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(sub_ren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(sub_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(smooth_sweeps_) = 1
p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_bjac_
p%baseprecv(ilev_)%iprcparm(mld_sub_solve_) = mld_ilu_n_
p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_ren_) = 0
p%baseprecv(ilev_)%iprcparm(mld_n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(mld_sub_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(mld_smooth_sweeps_) = 1
case ('AS')
nlev_ = 1
ilev_ = 1
allocate(p%baseprecv(nlev_),stat=info)
if (info == 0) call psb_realloc(ifpsz,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(dfpsz,p%baseprecv(ilev_)%dprcparm,info)
if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(mld_dfpsz_,p%baseprecv(ilev_)%dprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(prec_type_) = as_
p%baseprecv(ilev_)%iprcparm(sub_solve_) = ilu_n_
p%baseprecv(ilev_)%iprcparm(sub_restr_) = psb_halo_
p%baseprecv(ilev_)%iprcparm(sub_prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(sub_ren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(sub_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(smooth_sweeps_) = 1
p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_as_
p%baseprecv(ilev_)%iprcparm(mld_sub_solve_) = mld_ilu_n_
p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_halo_
p%baseprecv(ilev_)%iprcparm(mld_sub_prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_ren_) = 0
p%baseprecv(ilev_)%iprcparm(mld_n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(mld_sub_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(mld_smooth_sweeps_) = 1
case ('ML')
@ -136,61 +136,61 @@ subroutine mld_dprecinit(p,ptype,info,nlev)
endif
ilev_ = 1
allocate(p%baseprecv(nlev_),stat=info)
if (info == 0) call psb_realloc(ifpsz,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(dfpsz,p%baseprecv(ilev_)%dprcparm,info)
if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(mld_dfpsz_,p%baseprecv(ilev_)%dprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(prec_type_) = as_
p%baseprecv(ilev_)%iprcparm(sub_solve_) = ilu_n_
p%baseprecv(ilev_)%iprcparm(sub_restr_) = psb_halo_
p%baseprecv(ilev_)%iprcparm(sub_prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(sub_ren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(sub_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(smooth_sweeps_) = 1
p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_as_
p%baseprecv(ilev_)%iprcparm(mld_sub_solve_) = mld_ilu_n_
p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_halo_
p%baseprecv(ilev_)%iprcparm(mld_sub_prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_ren_) = 0
p%baseprecv(ilev_)%iprcparm(mld_n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(mld_sub_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(mld_smooth_sweeps_) = 1
if (nlev_ == 1) return
do ilev_ = 2, nlev_ -1
if (info == 0) call psb_realloc(ifpsz,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(dfpsz,p%baseprecv(ilev_)%dprcparm,info)
if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(mld_dfpsz_,p%baseprecv(ilev_)%dprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(prec_type_) = bjac_
p%baseprecv(ilev_)%iprcparm(sub_restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(sub_prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(sub_ren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(ml_type_) = mult_ml
p%baseprecv(ilev_)%iprcparm(aggr_alg_) = dec_aggr_
p%baseprecv(ilev_)%iprcparm(aggr_kind_) = smooth_prol_
p%baseprecv(ilev_)%iprcparm(coarse_mat_) = distr_mat_
p%baseprecv(ilev_)%iprcparm(smooth_pos_) = post_smooth_
p%baseprecv(ilev_)%iprcparm(aggr_eig_) = max_norm_
p%baseprecv(ilev_)%iprcparm(sub_solve_) = ilu_n_
p%baseprecv(ilev_)%iprcparm(sub_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(smooth_sweeps_) = 1
p%baseprecv(ilev_)%dprcparm(aggr_damp_) = 4.d0/3.d0
p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_bjac_
p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_ren_) = 0
p%baseprecv(ilev_)%iprcparm(mld_n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(mld_ml_type_) = mld_mult_ml_
p%baseprecv(ilev_)%iprcparm(mld_aggr_alg_) = mld_dec_aggr_
p%baseprecv(ilev_)%iprcparm(mld_aggr_kind_) = mld_smooth_prol_
p%baseprecv(ilev_)%iprcparm(mld_coarse_mat_) = mld_distr_mat_
p%baseprecv(ilev_)%iprcparm(mld_smooth_pos_) = mld_post_smooth_
p%baseprecv(ilev_)%iprcparm(mld_aggr_eig_) = mld_max_norm_
p%baseprecv(ilev_)%iprcparm(mld_sub_solve_) = mld_ilu_n_
p%baseprecv(ilev_)%iprcparm(mld_sub_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(mld_smooth_sweeps_) = 1
p%baseprecv(ilev_)%dprcparm(mld_aggr_damp_) = 4.d0/3.d0
end do
ilev_ = nlev_
if (info == 0) call psb_realloc(ifpsz,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(dfpsz,p%baseprecv(ilev_)%dprcparm,info)
if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(mld_dfpsz_,p%baseprecv(ilev_)%dprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(prec_type_) = bjac_
p%baseprecv(ilev_)%iprcparm(sub_restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(sub_prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(sub_ren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(ml_type_) = mult_ml
p%baseprecv(ilev_)%iprcparm(aggr_alg_) = dec_aggr_
p%baseprecv(ilev_)%iprcparm(aggr_kind_) = smooth_prol_
p%baseprecv(ilev_)%iprcparm(coarse_mat_) = distr_mat_
p%baseprecv(ilev_)%iprcparm(smooth_pos_) = post_smooth_
p%baseprecv(ilev_)%iprcparm(aggr_eig_) = max_norm_
p%baseprecv(ilev_)%iprcparm(sub_solve_) = umf_
p%baseprecv(ilev_)%iprcparm(sub_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(smooth_sweeps_) = 4
p%baseprecv(ilev_)%dprcparm(aggr_damp_) = 4.d0/3.d0
p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_bjac_
p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_ren_) = 0
p%baseprecv(ilev_)%iprcparm(mld_n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(mld_ml_type_) = mld_mult_ml_
p%baseprecv(ilev_)%iprcparm(mld_aggr_alg_) = mld_dec_aggr_
p%baseprecv(ilev_)%iprcparm(mld_aggr_kind_) = mld_smooth_prol_
p%baseprecv(ilev_)%iprcparm(mld_coarse_mat_) = mld_distr_mat_
p%baseprecv(ilev_)%iprcparm(mld_smooth_pos_) = mld_post_smooth_
p%baseprecv(ilev_)%iprcparm(mld_aggr_eig_) = mld_max_norm_
p%baseprecv(ilev_)%iprcparm(mld_sub_solve_) = mld_umf_
p%baseprecv(ilev_)%iprcparm(mld_sub_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(mld_smooth_sweeps_) = 4
p%baseprecv(ilev_)%dprcparm(mld_aggr_damp_) = 4.d0/3.d0
case default
write(0,*) 'Unknown preconditioner type request "',ptype,'"'

42
mlprec/mld_dprecset.f90
Unescape Escape View File

@ -79,7 +79,7 @@ subroutine mld_dprecseti(p,what,val,info,ilev)
if (ilev_ == 1) then
! Rules for fine level are slightly different.
select case(what)
case(prec_type_,sub_solve_,sub_restr_,sub_prol_,sub_ren_,n_ovr_,sub_fill_in_,smooth_sweeps_)
case(mld_prec_type_,mld_sub_solve_,mld_sub_restr_,mld_sub_prol_,mld_sub_ren_,mld_n_ovr_,mld_sub_fill_in_,mld_smooth_sweeps_)
p%baseprecv(ilev_)%iprcparm(what) = val
case default
write(0,*) 'Error: trying to call PRECSET with an invalid WHAT'
@ -87,31 +87,31 @@ subroutine mld_dprecseti(p,what,val,info,ilev)
end select
else if (ilev_ > 1) then
select case(what)
case(prec_type_,sub_solve_,sub_restr_,sub_prol_,sub_ren_,n_ovr_,sub_fill_in_,&
& smooth_sweeps_,ml_type_,aggr_alg_,aggr_kind_,coarse_mat_,&
& smooth_pos_,aggr_eig_)
case(mld_prec_type_,mld_sub_solve_,mld_sub_restr_,mld_sub_prol_,mld_sub_ren_,mld_n_ovr_,mld_sub_fill_in_,&
& mld_smooth_sweeps_,mld_ml_type_,mld_aggr_alg_,mld_aggr_kind_,mld_coarse_mat_,&
& mld_smooth_pos_,mld_aggr_eig_)
p%baseprecv(ilev_)%iprcparm(what) = val
case(coarse_solve_)
case(mld_coarse_solve_)
if (ilev_ /= nlev_) then
write(0,*) 'Inconsistent specification of WHAT vs. ILEV'
info = -2
return
end if
p%baseprecv(ilev_)%iprcparm(sub_solve_) = val
case(coarse_sweeps_)
p%baseprecv(ilev_)%iprcparm(mld_sub_solve_) = val
case(mld_coarse_sweeps_)
if (ilev_ /= nlev_) then
write(0,*) 'Inconsistent specification of WHAT vs. ILEV'
info = -2
return
end if
p%baseprecv(ilev_)%iprcparm(smooth_sweeps_) = val
case(coarse_fill_in_)
p%baseprecv(ilev_)%iprcparm(mld_smooth_sweeps_) = val
case(mld_coarse_fill_in_)
if (ilev_ /= nlev_) then
write(0,*) 'Inconsistent specification of WHAT vs. ILEV'
info = -2
return
end if
p%baseprecv(ilev_)%iprcparm(sub_fill_in_) = val
p%baseprecv(ilev_)%iprcparm(mld_sub_fill_in_) = val
case default
write(0,*) 'Error: trying to call PRECSET with an invalid WHAT'
info = -2
@ -121,9 +121,9 @@ subroutine mld_dprecseti(p,what,val,info,ilev)
else if (.not.present(ilev)) then
select case(what)
case(prec_type_,sub_solve_,sub_restr_,sub_prol_,sub_ren_,n_ovr_,sub_fill_in_,&
& smooth_sweeps_,ml_type_,aggr_alg_,aggr_kind_,coarse_mat_,&
& smooth_pos_,aggr_eig_)
case(mld_prec_type_,mld_sub_solve_,mld_sub_restr_,mld_sub_prol_,mld_sub_ren_,mld_n_ovr_,mld_sub_fill_in_,&
& mld_smooth_sweeps_,mld_ml_type_,mld_aggr_alg_,mld_aggr_kind_,mld_coarse_mat_,&
& mld_smooth_pos_,mld_aggr_eig_)
do ilev_=1,nlev_-1
if (.not.allocated(p%baseprecv(ilev_)%iprcparm)) then
write(0,*) 'Error: trying to call PRECSET on an uninitialized preconditioner'
@ -132,27 +132,27 @@ subroutine mld_dprecseti(p,what,val,info,ilev)
endif
p%baseprecv(ilev_)%iprcparm(what) = val
end do
case(coarse_solve_)
case(mld_coarse_solve_)
if (.not.allocated(p%baseprecv(nlev_)%iprcparm)) then
write(0,*) 'Error: trying to call PRECSET on an uninitialized preconditioner'
info = -1
return
endif
p%baseprecv(nlev_)%iprcparm(sub_solve_) = val
case(coarse_sweeps_)
p%baseprecv(nlev_)%iprcparm(mld_sub_solve_) = val
case(mld_coarse_sweeps_)
if (.not.allocated(p%baseprecv(nlev_)%iprcparm)) then
write(0,*) 'Error: trying to call PRECSET on an uninitialized preconditioner'
info = -1
return
endif
p%baseprecv(nlev_)%iprcparm(smooth_sweeps_) = val
case(coarse_fill_in_)
p%baseprecv(nlev_)%iprcparm(mld_smooth_sweeps_) = val
case(mld_coarse_fill_in_)
if (.not.allocated(p%baseprecv(nlev_)%iprcparm)) then
write(0,*) 'Error: trying to call PRECSET on an uninitialized preconditioner'
info = -1
return
endif
p%baseprecv(nlev_)%iprcparm(sub_fill_in_) = val
p%baseprecv(nlev_)%iprcparm(mld_sub_fill_in_) = val
case default
write(0,*) 'Error: trying to call PRECSET with an invalid WHAT'
info = -2
@ -205,14 +205,14 @@ subroutine mld_dprecsetd(p,what,val,info,ilev)
! Rules for fine level are slightly different.
select case(what)
! Right now we don't have any at base level. Will change when
! we implement ilu_t_
! we implement mld_ilu_t_
case default
write(0,*) 'Error: trying to call PRECSET with an invalid WHAT'
info = -2
end select
else if (ilev_ > 1) then
select case(what)
case(aggr_damp_)
case(mld_aggr_damp_)
p%baseprecv(ilev_)%dprcparm(what) = val
case default
write(0,*) 'Error: trying to call PRECSET with an invalid WHAT'

10
mlprec/mld_dslu_bld.f90
Unescape Escape View File

@ -59,7 +59,7 @@ subroutine mld_dslu_bld(a,desc_a,p,info)
call psb_info(ictxt, me, np)
if (toupper(a%fida) /= 'CSR') then
write(0,*) 'Unimplemented input to SLU_BLD'
write(0,*) 'Unimplemented input to mld_slu_BLD'
goto 9999
endif
@ -67,22 +67,22 @@ subroutine mld_dslu_bld(a,desc_a,p,info)
nzt = psb_sp_get_nnzeros(a)
if (Debug) then
write(0,*) me,'Calling psb_slu_factor ',nzt,a%m,&
write(0,*) me,'Calling psb_mld_slu_factor ',nzt,a%m,&
& a%k,p%desc_data%matrix_data(psb_n_row_)
call psb_barrier(ictxt)
endif
call mld_dslu_factor(a%m,nzt,&
& a%aspk,a%ia2,a%ia1,p%iprcparm(slu_ptr_),info)
& a%aspk,a%ia2,a%ia1,p%iprcparm(mld_slu_ptr_),info)
if (info /= 0) then
ch_err='psb_slu_fact'
ch_err='psb_mld_slu_fact'
call psb_errpush(4110,name,a_err=ch_err,i_err=(/info,0,0,0,0/))
goto 9999
end if
if (Debug) then
write(0,*) me, 'SPLUBLD: Done slu_Factor',info,p%iprcparm(slu_ptr_)
write(0,*) me, 'SPLUBLD: Done mld_slu_Factor',info,p%iprcparm(mld_slu_ptr_)
call psb_barrier(ictxt)
endif

4
mlprec/mld_dslud_bld.f90
Unescape Escape View File

@ -60,7 +60,7 @@ subroutine mld_dsludist_bld(a,desc_a,p,info)
call psb_info(ictxt, me, np)
if (toupper(a%fida) /= 'CSR') then
write(0,*) 'Unimplemented input to SLU_BLD'
write(0,*) 'Unimplemented input to mld_slu_BLD'
goto 9999
endif
@ -86,7 +86,7 @@ subroutine mld_dsludist_bld(a,desc_a,p,info)
call psb_loc_to_glob(a%ia1(1:nzt),desc_a,info,iact='I')
call mld_dsludist_factor(mglob,nrow,nzt,ifrst,&
& a%aspk,a%ia2,a%ia1,p%iprcparm(slud_ptr_),&
& a%aspk,a%ia2,a%ia1,p%iprcparm(mld_slud_ptr_),&
& npr, npc, info)
if (info /= 0) then
ch_err='psb_slud_fact'

6
mlprec/mld_dsp_renum.f90
Unescape Escape View File

@ -76,7 +76,7 @@ subroutine mld_dsp_renum(a,desc_a,blck,p,atmp,info)
call psb_spcnv(a,atmp,info,afmt='coo',dupl=psb_dupl_add_)
call psb_rwextd(a%m+blck%m,atmp,info,blck)
if (p%iprcparm(sub_ren_)==renum_glb_) then
if (p%iprcparm(mld_sub_ren_)==mld_renum_glb_) then
! This is the renumbering coherent with global indices..
mglob = psb_cd_get_global_rows(desc_a)
@ -106,7 +106,7 @@ subroutine mld_dsp_renum(a,desc_a,blck,p,atmp,info)
enddo
t3 = psb_wtime()
else if (p%iprcparm(sub_ren_)==renum_gps_) then
else if (p%iprcparm(mld_sub_ren_)==mld_renum_gps_) then
call psb_spcnv(atmp,info,afmt='csr',dupl=psb_dupl_add_)
nztmp = psb_sp_get_nnzeros(atmp)
@ -140,7 +140,7 @@ subroutine mld_dsp_renum(a,desc_a,blck,p,atmp,info)
itmp(1:8) = 0
! write(0,*) me,' Renumbering: Calling Metis'
! write(0,*) size(p%av(u_pr_)%pl),size(p%av(l_pr_)%pr)
! write(0,*) size(p%av(mld_u_pr_)%pl),size(p%av(mld_l_pr_)%pr)
call gps_reduction(atmp%m,atmp%ia2,atmp%ia1,p%perm,p%invperm,info)
if(info/=0) then
info=4010

10
mlprec/mld_dumf_bld.f90
Unescape Escape View File

@ -58,7 +58,7 @@ subroutine mld_dumf_bld(a,desc_a,p,info)
call psb_info(ictxt, me, np)
if (toupper(a%fida) /= 'CSC') then
write(0,*) 'Unimplemented input to UMF_BLD'
write(0,*) 'Unimplemented input to mld_umf_BLD'
goto 9999
endif
@ -66,7 +66,7 @@ subroutine mld_dumf_bld(a,desc_a,p,info)
nzt = psb_sp_get_nnzeros(a)
if (Debug) then
write(0,*) me,'Calling psb_umf_factor ',nzt,a%m,&
write(0,*) me,'Calling psb_mld_umf_factor ',nzt,a%m,&
& a%k,p%desc_data%matrix_data(psb_n_row_)
open(80+me)
call psb_csprt(80+me,a)
@ -76,17 +76,17 @@ subroutine mld_dumf_bld(a,desc_a,p,info)
call mld_dumf_factor(a%m,nzt,&
& a%aspk,a%ia1,a%ia2,&
& p%iprcparm(umf_symptr_),p%iprcparm(umf_numptr_),info)
& p%iprcparm(mld_umf_symptr_),p%iprcparm(mld_umf_numptr_),info)
if (info /= 0) then
i_err(1) = info
info=4110
call psb_errpush(info,name,a_err='psb_umf_fact',i_err=i_err)
call psb_errpush(info,name,a_err='psb_mld_umf_fact',i_err=i_err)
goto 9999
end if
if (Debug) then
write(0,*) me, 'UMFBLD: Done umf_Factor',info,p%iprcparm(umf_numptr_)
write(0,*) me, 'UMFBLD: Done mld_umf_Factor',info,p%iprcparm(mld_umf_numptr_)
call psb_barrier(ictxt)
endif

400
mlprec/mld_prec_type.f90
Unescape Escape View File

@ -65,11 +65,11 @@ module mld_prec_type
! 1. Number of levels = NLEV = size(baseprecv(:))
! 2. baseprecv(ilev)%av(:) sparse matrices needed for the current level.
! Includes:
! 2.1.: baseprecv(ilev)%av(l_pr_) L factor of ILU preconditioners
! 2.2.: baseprecv(ilev)%av(u_pr_) U factor of ILU preconditioners
! 2.3.: baseprecv(ilev)%av(ap_nd_) Off-diagonal part of A for Jacobi sweeps
! 2.4.: baseprecv(ilev)%av(ac_) Aggregated matrix of level ILEV
! 2.5.: baseprecv(ilev)%av(sm_pr_t_) Smoother prolongator transpose; maps vectors
! 2.1.: baseprecv(ilev)%av(mld_l_pr_) L factor of ILU preconditioners
! 2.2.: baseprecv(ilev)%av(mld_u_pr_) U factor of ILU preconditioners
! 2.3.: baseprecv(ilev)%av(mld_ap_nd_) Off-diagonal part of A for Jacobi sweeps
! 2.4.: baseprecv(ilev)%av(mld_ac_) Aggregated matrix of level ILEV
! 2.5.: baseprecv(ilev)%av(mld_sm_pr_t_) Smoother prolongator transpose; maps vectors
! (ilev-1) ---> (ilev)
! 2.6.: baseprecv(ilev)%av(sm_pr_) Smoother prolongator; maps vectors
! (ilev) ---> (ilev-1)
@ -78,7 +78,7 @@ module mld_prec_type
! 3. baseprecv(ilev)%desc_data comm descriptor for level ILEV
! 4. baseprecv(ilev)%base_a Pointer (really a pointer!) to the base matrix
! of the current level, i.e.: if ILEV=1 then A
! else the aggregated matrix av(ac_); so we have
! else the aggregated matrix av(mld_ac_); so we have
! a unified treatment of residuals. Need this to
! avoid passing explicitly matrix A to the
! outer prec. routine
@ -132,67 +132,67 @@ module mld_prec_type
! Entries in iprcparm
integer, parameter :: prec_type_=1
integer, parameter :: sub_solve_=2
integer, parameter :: sub_restr_=3
integer, parameter :: sub_prol_=4
integer, parameter :: sub_ren_=5
integer, parameter :: n_ovr_=6
integer, parameter :: sub_fill_in_=8
integer, parameter :: smooth_sweeps_=9
integer, parameter :: ml_type_=10
integer, parameter :: smooth_pos_=11
integer, parameter :: aggr_alg_=12
integer, parameter :: aggr_kind_=13
integer, parameter :: aggr_eig_=14
integer, parameter :: coarse_mat_=16
integer, parameter :: mld_prec_type_=1
integer, parameter :: mld_sub_solve_=2
integer, parameter :: mld_sub_restr_=3
integer, parameter :: mld_sub_prol_=4
integer, parameter :: mld_sub_ren_=5
integer, parameter :: mld_n_ovr_=6
integer, parameter :: mld_sub_fill_in_=8
integer, parameter :: mld_smooth_sweeps_=9
integer, parameter :: mld_ml_type_=10
integer, parameter :: mld_smooth_pos_=11
integer, parameter :: mld_aggr_alg_=12
integer, parameter :: mld_aggr_kind_=13
integer, parameter :: mld_aggr_eig_=14
integer, parameter :: mld_coarse_mat_=16
!! 2 ints for 64 bit versions
integer, parameter :: slu_ptr_=17
integer, parameter :: umf_symptr_=17
integer, parameter :: umf_numptr_=19
integer, parameter :: slud_ptr_=21
integer, parameter :: prec_status_=24
integer, parameter :: coarse_solve_ =25
integer, parameter :: coarse_sweeps_ =26
integer, parameter :: coarse_fill_in_=27
integer, parameter :: ifpsz=32
! Legal values for entry: prec_type_
integer, parameter :: min_prec_=0, noprec_=0, diag_=1, bjac_=2,&
& as_=3, max_prec_=3
! Legal values for entry: ml_type_
integer, parameter :: no_ml_=0, add_ml_=1, mult_ml=2
integer, parameter :: new_ml_prec_=3, max_ml_=new_ml_prec_
! Legal values for entry: smooth_pos_
integer, parameter :: pre_smooth_=1, post_smooth_=2, twoside_smooth_=3,&
& max_smooth_=twoside_smooth_
! Legal values for entry: sub_solve_
integer, parameter :: f_none_=0,ilu_n_=1,milu_n_=2, ilu_t_=3
integer, parameter :: slu_=4, umf_=5, sludist_=6
! Legal values for entry: aggr_alg_
integer, parameter :: dec_aggr_=0, sym_dec_aggr_=1
integer, parameter :: glb_aggr_=2, new_dec_aggr_=3
integer, parameter :: new_glb_aggr_=4, max_aggr_=new_glb_aggr_
! Legal values for entry: aggr_kind_
integer, parameter :: no_smooth_=0, smooth_prol_=1, biz_prol_=2
! Legal values for entry: aggr_eig_
integer, parameter :: max_norm_=0, user_choice_=999
! Legal values for entry: coarse_mat_
integer, parameter :: distr_mat_=0, repl_mat_=1
! Legal values for entry: prec_status_
integer, parameter :: prec_built=98765
! Legal values for entry: sub_ren_
integer, parameter :: renum_none_=0, renum_glb_=1, renum_gps_=2
integer, parameter :: mld_slu_ptr_=17
integer, parameter :: mld_umf_symptr_=17
integer, parameter :: mld_umf_numptr_=19
integer, parameter :: mld_slud_ptr_=21
integer, parameter :: mld_prec_status_=24
integer, parameter :: mld_coarse_solve_ =25
integer, parameter :: mld_coarse_sweeps_ =26
integer, parameter :: mld_coarse_fill_in_=27
integer, parameter :: mld_ifpsz_=32
! Legal values for entry: mld_prec_type_
integer, parameter :: mld_min_prec_=0, mld_noprec_=0, mld_diag_=1, mld_bjac_=2,&
& mld_as_=3, mld_max_prec_=3
! Legal values for entry: mld_ml_type_
integer, parameter :: mld_no_ml_=0, mld_add_ml_=1, mld_mult_ml_=2
integer, parameter :: mld_new_ml_prec_=3, mld_max_ml_=mld_new_ml_prec_
! Legal values for entry: mld_smooth_pos_
integer, parameter :: mld_pre_smooth_=1, mld_post_smooth_=2, mld_twoside_smooth_=3,&
& mld_max_smooth_=mld_twoside_smooth_
! Legal values for entry: mld_sub_solve_
integer, parameter :: mld_f_none_=0,mld_ilu_n_=1,mld_milu_n_=2, mld_ilu_t_=3
integer, parameter :: mld_slu_=4, mld_umf_=5, mld_sludist_=6
! Legal values for entry: mld_aggr_alg_
integer, parameter :: mld_dec_aggr_=0, sym_mld_dec_aggr_=1
integer, parameter :: mld_glb_aggr_=2, new_mld_dec_aggr_=3
integer, parameter :: new_mld_glb_aggr_=4, mld_max_aggr_=new_mld_glb_aggr_
! Legal values for entry: mld_aggr_kind_
integer, parameter :: mld_no_smooth_=0, mld_smooth_prol_=1, mld_biz_prol_=2
! Legal values for entry: mld_aggr_eig_
integer, parameter :: mld_max_norm_=0, mld_user_choice_=999
! Legal values for entry: mld_coarse_mat_
integer, parameter :: mld_distr_mat_=0, mld_repl_mat_=1
! Legal values for entry: mld_prec_status_
integer, parameter :: mld_prec_built_=98765
! Legal values for entry: mld_sub_ren_
integer, parameter :: mld_renum_none_=0, mld_renum_glb_=1, mld_renum_gps_=2
! Entries in dprcparm: ILU(T) epsilon, smoother omega
integer, parameter :: fact_eps_=1
integer, parameter :: aggr_damp_=2
integer, parameter :: aggr_thresh_=3
integer, parameter :: dfpsz=4
integer, parameter :: mld_fact_eps_=1
integer, parameter :: mld_aggr_damp_=2
integer, parameter :: mld_aggr_thresh_=3
integer, parameter :: mld_dfpsz_=4
! Fields for sparse matrices ensembles stored in av()
integer, parameter :: l_pr_=1, u_pr_=2, bp_ilu_avsz=2
integer, parameter :: ap_nd_=3, ac_=4, sm_pr_t_=5, sm_pr_=6
integer, parameter :: smth_avsz=6, max_avsz=smth_avsz
integer, parameter :: mld_l_pr_=1, mld_u_pr_=2, mld_bp_ilu_avsz_=2
integer, parameter :: mld_ap_nd_=3, mld_ac_=4, mld_sm_pr_t_=5, sm_pr_=6
integer, parameter :: mld_smth_avsz_=6, mld_max_avsz_=mld_smth_avsz_
@ -287,15 +287,15 @@ contains
val = 0
if (allocated(prec%iprcparm)) then
val = val + 4 * size(prec%iprcparm)
if (prec%iprcparm(prec_status_) == prec_built) then
select case(prec%iprcparm(sub_solve_))
case(ilu_n_,ilu_t_)
if (prec%iprcparm(mld_prec_status_) == mld_prec_built_) then
select case(prec%iprcparm(mld_sub_solve_))
case(mld_ilu_n_,mld_ilu_t_)
! do nothing
case(slu_)
case(mld_slu_)
write(0,*) 'Should implement check for size of SuperLU data structs'
case(umf_)
case(mld_umf_)
write(0,*) 'Should implement check for size of UMFPACK data structs'
case(sludist_)
case(mld_sludist_)
write(0,*) 'Should implement check for size of SuperLUDist data structs'
case default
end select
@ -326,15 +326,15 @@ contains
val = 0
if (allocated(prec%iprcparm)) then
val = val + 4 * size(prec%iprcparm)
if (prec%iprcparm(prec_status_) == prec_built) then
select case(prec%iprcparm(sub_solve_))
case(ilu_n_,ilu_t_)
if (prec%iprcparm(mld_prec_status_) == mld_prec_built_) then
select case(prec%iprcparm(mld_sub_solve_))
case(mld_ilu_n_,mld_ilu_t_)
! do nothing
case(slu_)
case(mld_slu_)
write(0,*) 'Should implement check for size of SuperLU data structs'
case(umf_)
case(mld_umf_)
write(0,*) 'Should implement check for size of UMFPACK data structs'
case(sludist_)
case(mld_sludist_)
write(0,*) 'Should implement check for size of SuperLUDist data structs'
case default
end select
@ -380,23 +380,23 @@ contains
if (allocated(p%baseprecv)) then
if (size(p%baseprecv)>=1) then
write(iout,*) 'Base preconditioner'
select case(p%baseprecv(1)%iprcparm(prec_type_))
case(noprec_)
select case(p%baseprecv(1)%iprcparm(mld_prec_type_))
case(mld_noprec_)
write(iout,*) 'No preconditioning'
case(diag_)
case(mld_diag_)
write(iout,*) 'Diagonal scaling'
case(bjac_)
case(mld_bjac_)
write(iout,*) 'Block Jacobi with: ',&
& fact_names(p%baseprecv(1)%iprcparm(sub_solve_))
case(as_)
& fact_names(p%baseprecv(1)%iprcparm(mld_sub_solve_))
case(mld_as_)
write(iout,*) 'Additive Schwarz with: ',&
& fact_names(p%baseprecv(1)%iprcparm(sub_solve_))
& fact_names(p%baseprecv(1)%iprcparm(mld_sub_solve_))
write(iout,*) 'Overlap:',&
& p%baseprecv(1)%iprcparm(n_ovr_)
& p%baseprecv(1)%iprcparm(mld_n_ovr_)
write(iout,*) 'Restriction: ',&
& restrict_names(p%baseprecv(1)%iprcparm(sub_restr_))
& restrict_names(p%baseprecv(1)%iprcparm(mld_sub_restr_))
write(iout,*) 'Prolongation: ',&
& prolong_names(p%baseprecv(1)%iprcparm(sub_prol_))
& prolong_names(p%baseprecv(1)%iprcparm(mld_sub_prol_))
end select
end if
if (size(p%baseprecv)>=2) then
@ -408,37 +408,37 @@ contains
write(iout,*) 'Multilevel: Level No', ilev
write(iout,*) 'Multilevel type: ',&
& ml_names(p%baseprecv(ilev)%iprcparm(ml_type_))
if (p%baseprecv(ilev)%iprcparm(ml_type_)>no_ml_) then
& ml_names(p%baseprecv(ilev)%iprcparm(mld_ml_type_))
if (p%baseprecv(ilev)%iprcparm(mld_ml_type_)>mld_no_ml_) then
write(iout,*) 'Multilevel aggregation: ', &
& aggr_names(p%baseprecv(ilev)%iprcparm(aggr_alg_))
& aggr_names(p%baseprecv(ilev)%iprcparm(mld_aggr_alg_))
write(iout,*) 'Smoother: ', &
& smooth_kinds(p%baseprecv(ilev)%iprcparm(aggr_kind_))
if (p%baseprecv(ilev)%iprcparm(aggr_kind_) /= no_smooth_) then
& smooth_kinds(p%baseprecv(ilev)%iprcparm(mld_aggr_kind_))
if (p%baseprecv(ilev)%iprcparm(mld_aggr_kind_) /= mld_no_smooth_) then
write(iout,*) 'Smoothing omega: ', &
& p%baseprecv(ilev)%dprcparm(aggr_damp_)
& p%baseprecv(ilev)%dprcparm(mld_aggr_damp_)
write(iout,*) 'Smoothing position: ',&
& smooth_names(p%baseprecv(ilev)%iprcparm(smooth_pos_))
& smooth_names(p%baseprecv(ilev)%iprcparm(mld_smooth_pos_))
end if
write(iout,*) 'Coarse matrix: ',&
& matrix_names(p%baseprecv(ilev)%iprcparm(coarse_mat_))
& matrix_names(p%baseprecv(ilev)%iprcparm(mld_coarse_mat_))
if (allocated(p%baseprecv(ilev)%nlaggr)) then
write(iout,*) 'Aggregation sizes: ', &
& sum( p%baseprecv(ilev)%nlaggr(:)),' : ',p%baseprecv(ilev)%nlaggr(:)
end if
write(iout,*) 'Factorization type: ',&
& fact_names(p%baseprecv(ilev)%iprcparm(sub_solve_))
select case(p%baseprecv(ilev)%iprcparm(sub_solve_))
case(ilu_n_)
write(iout,*) 'Fill level :',p%baseprecv(ilev)%iprcparm(sub_fill_in_)
case(ilu_t_)
write(iout,*) 'Fill threshold :',p%baseprecv(ilev)%dprcparm(fact_eps_)
case(slu_,umf_,sludist_)
& fact_names(p%baseprecv(ilev)%iprcparm(mld_sub_solve_))
select case(p%baseprecv(ilev)%iprcparm(mld_sub_solve_))
case(mld_ilu_n_)
write(iout,*) 'Fill level :',p%baseprecv(ilev)%iprcparm(mld_sub_fill_in_)
case(mld_ilu_t_)
write(iout,*) 'Fill threshold :',p%baseprecv(ilev)%dprcparm(mld_fact_eps_)
case(mld_slu_,mld_umf_,mld_sludist_)
case default
write(iout,*) 'Should never get here!'
end select
write(iout,*) 'Number of Jacobi sweeps: ', &
& (p%baseprecv(ilev)%iprcparm(smooth_sweeps_))
& (p%baseprecv(ilev)%iprcparm(mld_smooth_sweeps_))
end if
end do
end if
@ -459,23 +459,23 @@ contains
!!$ if (associated(p%baseprecv)) then
!!$ if (size(p%baseprecv)>=1) then
!!$ write(iout,*) 'Base preconditioner'
!!$ select case(p%baseprecv(1)%iprcparm(prec_type_))
!!$ case(noprec_)
!!$ select case(p%baseprecv(1)%iprcparm(mld_prec_type_))
!!$ case(mld_noprec_)
!!$ write(iout,*) 'No preconditioning'
!!$ case(diag_)
!!$ case(mld_diag_)
!!$ write(iout,*) 'Diagonal scaling'
!!$ case(bjac_)
!!$ case(mld_bjac_)
!!$ write(iout,*) 'Block Jacobi with: ',&
!!$ & fact_names(p%baseprecv(1)%iprcparm(sub_solve_))
!!$ case(as_,ras_,ash_,rash_)
!!$ & fact_names(p%baseprecv(1)%iprcparm(mld_sub_solve_))
!!$ case(mld_as_,rmld_as_,ash_,rash_)
!!$ write(iout,*) 'Additive Schwarz with: ',&
!!$ & fact_names(p%baseprecv(1)%iprcparm(sub_solve_))
!!$ & fact_names(p%baseprecv(1)%iprcparm(mld_sub_solve_))
!!$ write(iout,*) 'Overlap:',&
!!$ & p%baseprecv(1)%iprcparm(n_ovr_)
!!$ & p%baseprecv(1)%iprcparm(mld_n_ovr_)
!!$ write(iout,*) 'Restriction: ',&
!!$ & restrict_names(p%baseprecv(1)%iprcparm(sub_restr_))
!!$ & restrict_names(p%baseprecv(1)%iprcparm(mld_sub_restr_))
!!$ write(iout,*) 'Prolongation: ',&
!!$ & prolong_names(p%baseprecv(1)%iprcparm(sub_prol_))
!!$ & prolong_names(p%baseprecv(1)%iprcparm(mld_sub_prol_))
!!$ end select
!!$ end if
!!$ if (size(p%baseprecv)>=2) then
@ -483,30 +483,30 @@ contains
!!$ write(iout,*) 'Inconsistent MLPREC part!'
!!$ return
!!$ endif
!!$ write(iout,*) 'Multilevel: ',ml_names(p%baseprecv(2)%iprcparm(ml_type_))
!!$ if (p%baseprecv(2)%iprcparm(ml_type_)>no_ml_) then
!!$ write(iout,*) 'Multilevel: ',ml_names(p%baseprecv(2)%iprcparm(mld_ml_type_))
!!$ if (p%baseprecv(2)%iprcparm(mld_ml_type_)>mld_no_ml_) then
!!$ write(iout,*) 'Multilevel aggregation: ', &
!!$ & aggr_names(p%baseprecv(2)%iprcparm(aggr_alg_))
!!$ & aggr_names(p%baseprecv(2)%iprcparm(mld_aggr_alg_))
!!$ write(iout,*) 'Smoother: ', &
!!$ & smooth_kinds(p%baseprecv(2)%iprcparm(aggr_kind_))
!!$ write(iout,*) 'Smoothing omega: ', p%baseprecv(2)%dprcparm(aggr_damp_)
!!$ & smooth_kinds(p%baseprecv(2)%iprcparm(mld_aggr_kind_))
!!$ write(iout,*) 'Smoothing omega: ', p%baseprecv(2)%dprcparm(mld_aggr_damp_)
!!$ write(iout,*) 'Smoothing position: ',&
!!$ & smooth_names(p%baseprecv(2)%iprcparm(smooth_pos_))
!!$ & smooth_names(p%baseprecv(2)%iprcparm(mld_smooth_pos_))
!!$ write(iout,*) 'Coarse matrix: ',&
!!$ & matrix_names(p%baseprecv(2)%iprcparm(coarse_mat_))
!!$ & matrix_names(p%baseprecv(2)%iprcparm(mld_coarse_mat_))
!!$ write(iout,*) 'Factorization type: ',&
!!$ & fact_names(p%baseprecv(2)%iprcparm(sub_solve_))
!!$ select case(p%baseprecv(2)%iprcparm(sub_solve_))
!!$ case(ilu_n_)
!!$ write(iout,*) 'Fill level :',p%baseprecv(2)%iprcparm(sub_fill_in_)
!!$ case(ilu_t_)
!!$ write(iout,*) 'Fill threshold :',p%baseprecv(2)%dprcparm(fact_eps_)
!!$ case(slu_,umf_,sludist_)
!!$ & fact_names(p%baseprecv(2)%iprcparm(mld_sub_solve_))
!!$ select case(p%baseprecv(2)%iprcparm(mld_sub_solve_))
!!$ case(mld_ilu_n_)
!!$ write(iout,*) 'Fill level :',p%baseprecv(2)%iprcparm(mld_sub_fill_in_)
!!$ case(mld_ilu_t_)
!!$ write(iout,*) 'Fill threshold :',p%baseprecv(2)%dprcparm(mld_fact_eps_)
!!$ case(mld_slu_,mld_umf_,mld_sludist_)
!!$ case default
!!$ write(iout,*) 'Should never get here!'
!!$ end select
!!$ write(iout,*) 'Number of Jacobi sweeps: ', &
!!$ & (p%baseprecv(2)%iprcparm(smooth_sweeps_))
!!$ & (p%baseprecv(2)%iprcparm(mld_smooth_sweeps_))
!!$
!!$ end if
!!$ end if
@ -530,23 +530,23 @@ contains
if (allocated(p%baseprecv)) then
if (size(p%baseprecv)>=1) then
write(iout,*) 'Base preconditioner'
select case(p%baseprecv(1)%iprcparm(prec_type_))
case(noprec_)
select case(p%baseprecv(1)%iprcparm(mld_prec_type_))
case(mld_noprec_)
write(iout,*) 'No preconditioning'
case(diag_)
case(mld_diag_)
write(iout,*) 'Diagonal scaling'
case(bjac_)
case(mld_bjac_)
write(iout,*) 'Block Jacobi with: ',&
& fact_names(p%baseprecv(1)%iprcparm(sub_solve_))
case(as_)
& fact_names(p%baseprecv(1)%iprcparm(mld_sub_solve_))
case(mld_as_)
write(iout,*) 'Additive Schwarz with: ',&
& fact_names(p%baseprecv(1)%iprcparm(sub_solve_))
& fact_names(p%baseprecv(1)%iprcparm(mld_sub_solve_))
write(iout,*) 'Overlap:',&
& p%baseprecv(1)%iprcparm(n_ovr_)
& p%baseprecv(1)%iprcparm(mld_n_ovr_)
write(iout,*) 'Restriction: ',&
& restrict_names(p%baseprecv(1)%iprcparm(sub_restr_))
& restrict_names(p%baseprecv(1)%iprcparm(mld_sub_restr_))
write(iout,*) 'Prolongation: ',&
& prolong_names(p%baseprecv(1)%iprcparm(sub_prol_))
& prolong_names(p%baseprecv(1)%iprcparm(mld_sub_prol_))
end select
end if
if (size(p%baseprecv)>=2) then
@ -558,37 +558,37 @@ contains
write(iout,*) 'Multilevel: Level No', ilev
write(iout,*) 'Multilevel type: ',&
& ml_names(p%baseprecv(ilev)%iprcparm(ml_type_))
if (p%baseprecv(ilev)%iprcparm(ml_type_)>no_ml_) then
& ml_names(p%baseprecv(ilev)%iprcparm(mld_ml_type_))
if (p%baseprecv(ilev)%iprcparm(mld_ml_type_)>mld_no_ml_) then
write(iout,*) 'Multilevel aggregation: ', &
& aggr_names(p%baseprecv(ilev)%iprcparm(aggr_alg_))
& aggr_names(p%baseprecv(ilev)%iprcparm(mld_aggr_alg_))
write(iout,*) 'Smoother: ', &
& smooth_kinds(p%baseprecv(ilev)%iprcparm(aggr_kind_))
if (p%baseprecv(ilev)%iprcparm(aggr_kind_) /= no_smooth_) then
& smooth_kinds(p%baseprecv(ilev)%iprcparm(mld_aggr_kind_))
if (p%baseprecv(ilev)%iprcparm(mld_aggr_kind_) /= mld_no_smooth_) then
write(iout,*) 'Smoothing omega: ', &
& p%baseprecv(ilev)%dprcparm(aggr_damp_)
& p%baseprecv(ilev)%dprcparm(mld_aggr_damp_)
write(iout,*) 'Smoothing position: ',&
& smooth_names(p%baseprecv(ilev)%iprcparm(smooth_pos_))
& smooth_names(p%baseprecv(ilev)%iprcparm(mld_smooth_pos_))
end if
write(iout,*) 'Coarse matrix: ',&
& matrix_names(p%baseprecv(ilev)%iprcparm(coarse_mat_))
& matrix_names(p%baseprecv(ilev)%iprcparm(mld_coarse_mat_))
if (allocated(p%baseprecv(ilev)%nlaggr)) then
write(iout,*) 'Aggregation sizes: ', &
& sum( p%baseprecv(ilev)%nlaggr(:)),' : ',p%baseprecv(ilev)%nlaggr(:)
end if
write(iout,*) 'Factorization type: ',&
& fact_names(p%baseprecv(ilev)%iprcparm(sub_solve_))
select case(p%baseprecv(ilev)%iprcparm(sub_solve_))
case(ilu_n_)
write(iout,*) 'Fill level :',p%baseprecv(ilev)%iprcparm(sub_fill_in_)
case(ilu_t_)
write(iout,*) 'Fill threshold :',p%baseprecv(ilev)%dprcparm(fact_eps_)
case(slu_,umf_,sludist_)
& fact_names(p%baseprecv(ilev)%iprcparm(mld_sub_solve_))
select case(p%baseprecv(ilev)%iprcparm(mld_sub_solve_))
case(mld_ilu_n_)
write(iout,*) 'Fill level :',p%baseprecv(ilev)%iprcparm(mld_sub_fill_in_)
case(mld_ilu_t_)
write(iout,*) 'Fill threshold :',p%baseprecv(ilev)%dprcparm(mld_fact_eps_)
case(mld_slu_,mld_umf_,mld_sludist_)
case default
write(iout,*) 'Should never get here!'
end select
write(iout,*) 'Number of Jacobi sweeps: ', &
& (p%baseprecv(ilev)%iprcparm(smooth_sweeps_))
& (p%baseprecv(ilev)%iprcparm(mld_smooth_sweeps_))
end if
end do
end if
@ -609,23 +609,23 @@ contains
!!$ if (associated(p%baseprecv)) then
!!$ if (size(p%baseprecv)>=1) then
!!$ write(iout,*) 'Base preconditioner'
!!$ select case(p%baseprecv(1)%iprcparm(prec_type_))
!!$ case(noprec_)
!!$ select case(p%baseprecv(1)%iprcparm(mld_prec_type_))
!!$ case(mld_noprec_)
!!$ write(iout,*) 'No preconditioning'
!!$ case(diag_)
!!$ case(mld_diag_)
!!$ write(iout,*) 'Diagonal scaling'
!!$ case(bjac_)
!!$ case(mld_bjac_)
!!$ write(iout,*) 'Block Jacobi with: ',&
!!$ & fact_names(p%baseprecv(1)%iprcparm(sub_solve_))
!!$ case(as_,ras_,ash_,rash_)
!!$ & fact_names(p%baseprecv(1)%iprcparm(mld_sub_solve_))
!!$ case(mld_as_,rmld_as_,ash_,rash_)
!!$ write(iout,*) 'Additive Schwarz with: ',&
!!$ & fact_names(p%baseprecv(1)%iprcparm(sub_solve_))
!!$ & fact_names(p%baseprecv(1)%iprcparm(mld_sub_solve_))
!!$ write(iout,*) 'Overlap:',&
!!$ & p%baseprecv(1)%iprcparm(n_ovr_)
!!$ & p%baseprecv(1)%iprcparm(mld_n_ovr_)
!!$ write(iout,*) 'Restriction: ',&
!!$ & restrict_names(p%baseprecv(1)%iprcparm(sub_restr_))
!!$ & restrict_names(p%baseprecv(1)%iprcparm(mld_sub_restr_))
!!$ write(iout,*) 'Prolongation: ',&
!!$ & prolong_names(p%baseprecv(1)%iprcparm(sub_prol_))
!!$ & prolong_names(p%baseprecv(1)%iprcparm(mld_sub_prol_))
!!$ end select
!!$ end if
!!$ if (size(p%baseprecv)>=2) then
@ -633,30 +633,30 @@ contains
!!$ write(iout,*) 'Inconsistent MLPREC part!'
!!$ return
!!$ endif
!!$ write(iout,*) 'Multilevel: ',ml_names(p%baseprecv(2)%iprcparm(ml_type_))
!!$ if (p%baseprecv(2)%iprcparm(ml_type_)>no_ml_) then
!!$ write(iout,*) 'Multilevel: ',ml_names(p%baseprecv(2)%iprcparm(mld_ml_type_))
!!$ if (p%baseprecv(2)%iprcparm(mld_ml_type_)>mld_no_ml_) then
!!$ write(iout,*) 'Multilevel aggregation: ', &
!!$ & aggr_names(p%baseprecv(2)%iprcparm(aggr_alg_))
!!$ & aggr_names(p%baseprecv(2)%iprcparm(mld_aggr_alg_))
!!$ write(iout,*) 'Smoother: ', &
!!$ & smooth_kinds(p%baseprecv(2)%iprcparm(aggr_kind_))
!!$ write(iout,*) 'Smoothing omega: ', p%baseprecv(2)%dprcparm(aggr_damp_)
!!$ & smooth_kinds(p%baseprecv(2)%iprcparm(mld_aggr_kind_))
!!$ write(iout,*) 'Smoothing omega: ', p%baseprecv(2)%dprcparm(mld_aggr_damp_)
!!$ write(iout,*) 'Smoothing position: ',&
!!$ & smooth_names(p%baseprecv(2)%iprcparm(smooth_pos_))
!!$ & smooth_names(p%baseprecv(2)%iprcparm(mld_smooth_pos_))
!!$ write(iout,*) 'Coarse matrix: ',&
!!$ & matrix_names(p%baseprecv(2)%iprcparm(coarse_mat_))
!!$ & matrix_names(p%baseprecv(2)%iprcparm(mld_coarse_mat_))
!!$ write(iout,*) 'Factorization type: ',&
!!$ & fact_names(p%baseprecv(2)%iprcparm(sub_solve_))
!!$ select case(p%baseprecv(2)%iprcparm(sub_solve_))
!!$ case(ilu_n_)
!!$ write(iout,*) 'Fill level :',p%baseprecv(2)%iprcparm(sub_fill_in_)
!!$ case(ilu_t_)
!!$ write(iout,*) 'Fill threshold :',p%baseprecv(2)%dprcparm(fact_eps_)
!!$ case(slu_,umf_,sludist_)
!!$ & fact_names(p%baseprecv(2)%iprcparm(mld_sub_solve_))
!!$ select case(p%baseprecv(2)%iprcparm(mld_sub_solve_))
!!$ case(mld_ilu_n_)
!!$ write(iout,*) 'Fill level :',p%baseprecv(2)%iprcparm(mld_sub_fill_in_)
!!$ case(mld_ilu_t_)
!!$ write(iout,*) 'Fill threshold :',p%baseprecv(2)%dprcparm(mld_fact_eps_)
!!$ case(mld_slu_,mld_umf_,mld_sludist_)
!!$ case default
!!$ write(iout,*) 'Should never get here!'
!!$ end select
!!$ write(iout,*) 'Number of Jacobi sweeps: ', &
!!$ & (p%baseprecv(2)%iprcparm(smooth_sweeps_))
!!$ & (p%baseprecv(2)%iprcparm(mld_smooth_sweeps_))
!!$
!!$ end if
!!$ end if
@ -676,7 +676,7 @@ contains
integer, intent(in) :: ip
logical :: is_legal_base_prec
is_legal_base_prec = ((ip>=noprec_).and.(ip<=max_prec_))
is_legal_base_prec = ((ip>=mld_noprec_).and.(ip<=mld_max_prec_))
return
end function is_legal_base_prec
function is_legal_n_ovr(ip)
@ -722,7 +722,7 @@ contains
integer, intent(in) :: ip
logical :: is_legal_ml_type
is_legal_ml_type = ((ip>=no_ml_).and.(ip<=max_ml_))
is_legal_ml_type = ((ip>=mld_no_ml_).and.(ip<=mld_max_ml_))
return
end function is_legal_ml_type
function is_legal_ml_aggr_kind(ip)
@ -730,7 +730,7 @@ contains
integer, intent(in) :: ip
logical :: is_legal_ml_aggr_kind
is_legal_ml_aggr_kind = ((ip>=dec_aggr_).and.(ip<=max_aggr_))
is_legal_ml_aggr_kind = ((ip>=mld_dec_aggr_).and.(ip<=mld_max_aggr_))
return
end function is_legal_ml_aggr_kind
function is_legal_ml_smooth_pos(ip)
@ -738,7 +738,7 @@ contains
integer, intent(in) :: ip
logical :: is_legal_ml_smooth_pos
is_legal_ml_smooth_pos = ((ip>=pre_smooth_).and.(ip<=max_smooth_))
is_legal_ml_smooth_pos = ((ip>=mld_pre_smooth_).and.(ip<=mld_max_smooth_))
return
end function is_legal_ml_smooth_pos
function is_legal_ml_smth_kind(ip)
@ -746,7 +746,7 @@ contains
integer, intent(in) :: ip
logical :: is_legal_ml_smth_kind
is_legal_ml_smth_kind = ((ip>=no_smooth_).and.(ip<=biz_prol_))
is_legal_ml_smth_kind = ((ip>=mld_no_smooth_).and.(ip<=mld_biz_prol_))
return
end function is_legal_ml_smth_kind
function is_legal_ml_coarse_mat(ip)
@ -754,7 +754,7 @@ contains
integer, intent(in) :: ip
logical :: is_legal_ml_coarse_mat
is_legal_ml_coarse_mat = ((ip>=distr_mat_).and.(ip<=repl_mat_))
is_legal_ml_coarse_mat = ((ip>=mld_distr_mat_).and.(ip<=mld_repl_mat_))
return
end function is_legal_ml_coarse_mat
function is_legal_ml_fact(ip)
@ -762,7 +762,7 @@ contains
integer, intent(in) :: ip
logical :: is_legal_ml_fact
is_legal_ml_fact = ((ip>=ilu_n_).and.(ip<=sludist_))
is_legal_ml_fact = ((ip>=mld_ilu_n_).and.(ip<=mld_sludist_))
return
end function is_legal_ml_fact
function is_legal_ml_lev(ip)
@ -889,15 +889,15 @@ contains
endif
if (allocated(p%iprcparm)) then
if (p%iprcparm(sub_solve_)==slu_) then
call mld_dslu_free(p%iprcparm(slu_ptr_),info)
if (p%iprcparm(mld_sub_solve_)==mld_slu_) then
call mld_dslu_free(p%iprcparm(mld_slu_ptr_),info)
end if
if (p%iprcparm(sub_solve_)==sludist_) then
call mld_dsludist_free(p%iprcparm(slud_ptr_),info)
if (p%iprcparm(mld_sub_solve_)==mld_sludist_) then
call mld_dsludist_free(p%iprcparm(mld_slud_ptr_),info)
end if
if (p%iprcparm(sub_solve_)==umf_) then
call mld_dumf_free(p%iprcparm(umf_symptr_),&
& p%iprcparm(umf_numptr_),info)
if (p%iprcparm(mld_sub_solve_)==mld_umf_) then
call mld_dumf_free(p%iprcparm(mld_umf_symptr_),&
& p%iprcparm(mld_umf_numptr_),info)
end if
deallocate(p%iprcparm,stat=info)
end if
@ -972,12 +972,12 @@ contains
endif
if (allocated(p%iprcparm)) then
if (p%iprcparm(sub_solve_)==slu_) then
call mld_zslu_free(p%iprcparm(slu_ptr_),info)
if (p%iprcparm(mld_sub_solve_)==mld_slu_) then
call mld_zslu_free(p%iprcparm(mld_slu_ptr_),info)
end if
if (p%iprcparm(sub_solve_)==umf_) then
call mld_zumf_free(p%iprcparm(umf_symptr_),&
& p%iprcparm(umf_numptr_),info)
if (p%iprcparm(mld_sub_solve_)==mld_umf_) then
call mld_zumf_free(p%iprcparm(mld_umf_symptr_),&
& p%iprcparm(mld_umf_numptr_),info)
end if
deallocate(p%iprcparm,stat=info)
end if
@ -1003,13 +1003,13 @@ contains
character(len=10) :: pr_to_str
select case(iprec)
case(noprec_)
case(mld_noprec_)
pr_to_str='NOPREC'
case(diag_)
case(mld_diag_)
pr_to_str='DIAG'
case(bjac_)
case(mld_bjac_)
pr_to_str='BJAC'
case(as_)
case(mld_as_)
pr_to_str='AS'
end select

6
mlprec/mld_zaggrmap_bld.f90
Unescape Escape View File

@ -72,7 +72,7 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
ncol = psb_cd_get_local_cols(desc_a)
select case (aggr_type)
case (dec_aggr_,sym_dec_aggr_)
case (mld_dec_aggr_,sym_mld_dec_aggr_)
nr = a%m
@ -87,7 +87,7 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
do i=1, nr
ilaggr(i) = -(nr+1)
end do
if (aggr_type == dec_aggr_) then
if (aggr_type == mld_dec_aggr_) then
apnt => a
else
call psb_sp_clip(a,atmp,info,imax=nr,jmax=nr,&
@ -305,7 +305,7 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info)
nlaggr(me+1) = naggr
call psb_sum(ictxt,nlaggr(1:np))
if (aggr_type == sym_dec_aggr_) then
if (aggr_type == sym_mld_dec_aggr_) then
call psb_sp_free(atmp,info)
end if

6
mlprec/mld_zaggrmat_asb.f90
Unescape Escape View File

@ -61,8 +61,8 @@ subroutine mld_zaggrmat_asb(a,desc_a,ac,desc_ac,p,info)
call psb_info(ictxt, me, np)
select case (p%iprcparm(aggr_kind_))
case (no_smooth_)
select case (p%iprcparm(mld_aggr_kind_))
case (mld_no_smooth_)
call mld_aggrmat_raw_asb(a,desc_a,ac,desc_ac,p,info)
@ -72,7 +72,7 @@ subroutine mld_zaggrmat_asb(a,desc_a,ac,desc_ac,p,info)
end if
if (aggr_dump) call psb_csprt(90+me,ac,head='% Raw aggregate.')
case(smooth_prol_,biz_prol_)
case(mld_smooth_prol_,mld_biz_prol_)
if (aggr_dump) call psb_csprt(70+me,a,head='% Input matrix')
call mld_aggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)

8
mlprec/mld_zaggrmat_raw_asb.F90
Unescape Escape View File

@ -88,7 +88,7 @@ subroutine mld_zaggrmat_raw_asb(a,desc_a,ac,desc_ac,p,info)
naggrm1=sum(p%nlaggr(1:me))
if (p%iprcparm(coarse_mat_) == repl_mat_) then
if (p%iprcparm(mld_coarse_mat_) == mld_repl_mat_) then
do i=1, nrow
p%mlia(i) = p%mlia(i) + naggrm1
end do
@ -121,7 +121,7 @@ subroutine mld_zaggrmat_raw_asb(a,desc_a,ac,desc_ac,p,info)
! This is to minimize data exchange
call psb_spcnv(b,info,afmt='coo',dupl=psb_dupl_add_)
if (p%iprcparm(coarse_mat_) == repl_mat_) then
if (p%iprcparm(mld_coarse_mat_) == mld_repl_mat_) then
call psb_cdrep(ntaggr,ictxt,desc_ac,info)
if(info /= 0) then
@ -168,7 +168,7 @@ subroutine mld_zaggrmat_raw_asb(a,desc_a,ac,desc_ac,p,info)
goto 9999
end if
else if (p%iprcparm(coarse_mat_) == distr_mat_) then
else if (p%iprcparm(mld_coarse_mat_) == mld_distr_mat_) then
call psb_cdall(ictxt,desc_ac,info,nl=naggr)
if(info /= 0) then
@ -194,7 +194,7 @@ subroutine mld_zaggrmat_raw_asb(a,desc_a,ac,desc_ac,p,info)
else
write(0,*) 'Unknown p%iprcparm(coarse_mat) in aggregate_sp',p%iprcparm(coarse_mat_)
write(0,*) 'Unknown p%iprcparm(coarse_mat) in aggregate_sp',p%iprcparm(mld_coarse_mat_)
end if
deallocate(nzbr,idisp)

46
mlprec/mld_zaggrmat_smth_asb.F90
Unescape Escape View File

@ -84,7 +84,7 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
call psb_nullify_sp(am3)
call psb_nullify_sp(am4)
am2 => p%av(sm_pr_t_)
am2 => p%av(mld_sm_pr_t_)
am1 => p%av(sm_pr_)
call psb_nullify_sp(am1)
call psb_nullify_sp(am2)
@ -108,9 +108,9 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
naggrm1 = sum(p%nlaggr(1:me))
naggrp1 = sum(p%nlaggr(1:me+1))
ml_global_nmb = ( (p%iprcparm(aggr_kind_) == smooth_prol_).or.&
& ( (p%iprcparm(aggr_kind_) == biz_prol_).and.&
& (p%iprcparm(coarse_mat_) == repl_mat_)) )
ml_global_nmb = ( (p%iprcparm(mld_aggr_kind_) == mld_smooth_prol_).or.&
& ( (p%iprcparm(mld_aggr_kind_) == mld_biz_prol_).and.&
& (p%iprcparm(mld_coarse_mat_) == mld_repl_mat_)) )
if (ml_global_nmb) then
@ -217,9 +217,9 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
call psb_sp_scal(am3,p%dorig,info)
if(info /= 0) goto 9999
if (p%iprcparm(aggr_eig_) == max_norm_) then
if (p%iprcparm(mld_aggr_eig_) == mld_max_norm_) then
if (p%iprcparm(aggr_kind_) == biz_prol_) then
if (p%iprcparm(mld_aggr_kind_) == mld_biz_prol_) then
!
! This only works with CSR.
@ -244,15 +244,15 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
anorm = psb_spnrmi(am3,desc_a,info)
endif
omega = 4.d0/(3.d0*anorm)
p%dprcparm(aggr_damp_) = omega
p%dprcparm(mld_aggr_damp_) = omega
else if (p%iprcparm(aggr_eig_) == user_choice_) then
else if (p%iprcparm(mld_aggr_eig_) == mld_user_choice_) then
omega = p%dprcparm(aggr_damp_)
omega = p%dprcparm(mld_aggr_damp_)
else if (p%iprcparm(aggr_eig_) /= user_choice_) then
else if (p%iprcparm(mld_aggr_eig_) /= mld_user_choice_) then
write(0,*) me,'Error: invalid choice for OMEGA in blaggrmat?? ',&
& p%iprcparm(aggr_eig_)
& p%iprcparm(mld_aggr_eig_)
end if
@ -370,7 +370,7 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
call psb_numbmm(a,am1,am3)
if (debug) write(0,*) me,'Done NUMBMM 2'
if (p%iprcparm(aggr_kind_) == smooth_prol_) then
if (p%iprcparm(mld_aggr_kind_) == mld_smooth_prol_) then
call psb_transp(am1,am2,fmt='COO')
nzl = am2%infoa(psb_nnz_)
i=0
@ -398,7 +398,7 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
endif
if (debug) write(0,*) me,'starting sphalo/ rwxtd'
if (p%iprcparm(aggr_kind_) == smooth_prol_) then
if (p%iprcparm(mld_aggr_kind_) == mld_smooth_prol_) then
! am2 = ((i-wDA)Ptilde)^T
call psb_sphalo(am3,desc_a,am4,info,&
& colcnv=.false.,rowscale=.true.)
@ -418,7 +418,7 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
goto 9999
end if
else if (p%iprcparm(aggr_kind_) == biz_prol_) then
else if (p%iprcparm(mld_aggr_kind_) == mld_biz_prol_) then
call psb_rwextd(ncol,am3,info)
if(info /= 0) then
@ -454,13 +454,13 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
if (test_dump) call psb_csprt(80+me,b,head='% Smoothed aggregate AC.')
select case(p%iprcparm(aggr_kind_))
select case(p%iprcparm(mld_aggr_kind_))
case(smooth_prol_)
case(mld_smooth_prol_)
select case(p%iprcparm(coarse_mat_))
select case(p%iprcparm(mld_coarse_mat_))
case(distr_mat_)
case(mld_distr_mat_)
call psb_sp_clone(b,ac,info)
if(info /= 0) goto 9999
@ -550,7 +550,7 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
am2%m=desc_ac%matrix_data(psb_n_col_)
if (debug) write(0,*) me,'Done ac '
case(repl_mat_)
case(mld_repl_mat_)
!
!
call psb_cdrep(ntaggr,ictxt,desc_ac,info)
@ -602,11 +602,11 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
end select
case(biz_prol_)
case(mld_biz_prol_)
select case(p%iprcparm(coarse_mat_))
select case(p%iprcparm(mld_coarse_mat_))
case(distr_mat_)
case(mld_distr_mat_)
call psb_sp_clone(b,ac,info)
if(info /= 0) then
@ -632,7 +632,7 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ac,desc_ac,p,info)
end if
case(repl_mat_)
case(mld_repl_mat_)
!
!

4
mlprec/mld_zasmat_bld.f90
Unescape Escape View File

@ -98,7 +98,7 @@ Subroutine mld_zasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt)
nhalo = n_col-nrow_a
If (ptype == bjac_) Then
If (ptype == mld_bjac_) Then
!
! Block Jacobi. Copy the descriptor, just in case we want to
! do the renumbering.
@ -125,7 +125,7 @@ Subroutine mld_zasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt)
end if
endif
Else If (ptype == as_) Then
Else If (ptype == mld_as_) Then
!

32
mlprec/mld_zbaseprec_aply.f90
Unescape Escape View File

@ -82,13 +82,13 @@ subroutine mld_zbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
goto 9999
end select
select case(prec%iprcparm(prec_type_))
select case(prec%iprcparm(mld_prec_type_))
case(noprec_)
case(mld_noprec_)
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
case(diag_)
case(mld_diag_)
if (size(work) >= size(x)) then
ww => work
@ -112,7 +112,7 @@ subroutine mld_zbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
end if
end if
case(bjac_)
case(mld_bjac_)
call mld_bjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
if(info.ne.0) then
@ -121,9 +121,9 @@ subroutine mld_zbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
goto 9999
end if
case(as_)
case(mld_as_)
if (prec%iprcparm(n_ovr_)==0) then
if (prec%iprcparm(mld_n_ovr_)==0) then
! shortcut: this fixes performance for RAS(0) == BJA
call mld_bjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
if(info.ne.0) then
@ -176,19 +176,19 @@ subroutine mld_zbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
tx(1:nrow_d) = x(1:nrow_d)
tx(nrow_d+1:isz) = zzero
if (prec%iprcparm(sub_restr_)==psb_halo_) then
if (prec%iprcparm(mld_sub_restr_)==psb_halo_) then
call psb_halo(tx,prec%desc_data,info,work=aux,data=psb_comm_ext_)
if(info /=0) then
info=4010
ch_err='psb_halo'
goto 9999
end if
else if (prec%iprcparm(sub_restr_) /= psb_none_) then
else if (prec%iprcparm(mld_sub_restr_) /= psb_none_) then
write(0,*) 'Problem in PRC_APLY: Unknown value for restriction ',&
&prec%iprcparm(sub_restr_)
&prec%iprcparm(mld_sub_restr_)
end if
if (prec%iprcparm(sub_ren_)>0) then
if (prec%iprcparm(mld_sub_ren_)>0) then
call psb_gelp('n',prec%perm,tx,info)
!!$ call zgelp('N',n_row,1,prec%perm,tx,isz,ww,isz,info)
if(info /=0) then
@ -205,7 +205,7 @@ subroutine mld_zbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
goto 9999
end if
if (prec%iprcparm(sub_ren_)>0) then
if (prec%iprcparm(mld_sub_ren_)>0) then
call psb_gelp('n',prec%invperm,ty,info)
!!$ call zgelp('N',n_row,1,prec%invperm,ty,isz,ww,isz,info)
if(info /=0) then
@ -215,7 +215,7 @@ subroutine mld_zbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
end if
endif
select case (prec%iprcparm(sub_prol_))
select case (prec%iprcparm(mld_sub_prol_))
case(psb_none_)
! Would work anyway, but since it's supposed to do nothing...
@ -223,7 +223,7 @@ subroutine mld_zbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
case(psb_sum_,psb_avg_)
call psb_ovrl(ty,prec%desc_data,info,&
& update=prec%iprcparm(sub_prol_),work=aux)
& update=prec%iprcparm(mld_sub_prol_),work=aux)
if(info /=0) then
info=4010
ch_err='psb_ovrl'
@ -232,7 +232,7 @@ subroutine mld_zbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
case default
write(0,*) 'Problem in PRC_APLY: Unknown value for prolongation ',&
& prec%iprcparm(sub_prol_)
& prec%iprcparm(mld_sub_prol_)
end select
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
@ -248,8 +248,8 @@ subroutine mld_zbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
endif
end if
case default
write(0,*) 'Invalid PRE%PREC ',prec%iprcparm(prec_type_),':',&
& min_prec_,noprec_,diag_,bjac_,as_
write(0,*) 'Invalid PRE%PREC ',prec%iprcparm(mld_prec_type_),':',&
& mld_min_prec_,mld_noprec_,mld_diag_,mld_bjac_,mld_as_
end select
call psb_erractionrestore(err_act)

36
mlprec/mld_zbaseprec_bld.f90
Unescape Escape View File

@ -91,14 +91,14 @@ subroutine mld_zbaseprc_bld(a,desc_a,p,info,upd)
! Should add check to ensure all procs have the same...
!
call mld_check_def(p%iprcparm(prec_type_),'base_prec',&
& diag_,is_legal_base_prec)
call mld_check_def(p%iprcparm(mld_prec_type_),'base_prec',&
& mld_diag_,is_legal_base_prec)
call psb_nullify_desc(p%desc_data)
select case(p%iprcparm(prec_type_))
case (noprec_)
select case(p%iprcparm(mld_prec_type_))
case (mld_noprec_)
! Do nothing.
call psb_cdcpy(desc_a,p%desc_data,info)
if(info /= 0) then
@ -108,7 +108,7 @@ subroutine mld_zbaseprc_bld(a,desc_a,p,info,upd)
goto 9999
end if
case (diag_)
case (mld_diag_)
call mld_diag_bld(a,desc_a,p,iupd,info)
if(debug) write(0,*)me,': out of mld_diag_bld'
@ -119,22 +119,22 @@ subroutine mld_zbaseprc_bld(a,desc_a,p,info,upd)
goto 9999
end if
case (bjac_,as_)
case (mld_bjac_,mld_as_)
call mld_check_def(p%iprcparm(n_ovr_),'overlap',&
call mld_check_def(p%iprcparm(mld_n_ovr_),'overlap',&
& 0,is_legal_n_ovr)
call mld_check_def(p%iprcparm(sub_restr_),'restriction',&
call mld_check_def(p%iprcparm(mld_sub_restr_),'restriction',&
& psb_halo_,is_legal_restrict)
call mld_check_def(p%iprcparm(sub_prol_),'prolongator',&
call mld_check_def(p%iprcparm(mld_sub_prol_),'prolongator',&
& psb_none_,is_legal_prolong)
call mld_check_def(p%iprcparm(sub_ren_),'renumbering',&
& renum_none_,is_legal_renum)
call mld_check_def(p%iprcparm(sub_solve_),'fact',&
& ilu_n_,is_legal_ml_fact)
call mld_check_def(p%iprcparm(mld_sub_ren_),'renumbering',&
& mld_renum_none_,is_legal_renum)
call mld_check_def(p%iprcparm(mld_sub_solve_),'fact',&
& mld_ilu_n_,is_legal_ml_fact)
if (p%iprcparm(sub_solve_)==sludist_) then
p%iprcparm(n_ovr_) = 0
p%iprcparm(smooth_sweeps_) = 1
if (p%iprcparm(mld_sub_solve_)==mld_sludist_) then
p%iprcparm(mld_n_ovr_) = 0
p%iprcparm(mld_smooth_sweeps_) = 1
end if
if (debug) write(0,*)me, ': Calling mld_bjac_bld'
@ -149,7 +149,7 @@ subroutine mld_zbaseprc_bld(a,desc_a,p,info,upd)
case default
info=4010
ch_err='Unknown prec_type_'
ch_err='Unknown mld_prec_type_'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
@ -157,7 +157,7 @@ subroutine mld_zbaseprc_bld(a,desc_a,p,info,upd)
p%base_a => a
p%base_desc => desc_a
p%iprcparm(prec_status_) = prec_built
p%iprcparm(mld_prec_status_) = mld_prec_built_
call psb_erractionrestore(err_act)
return

80
mlprec/mld_zbjac_aply.f90
Unescape Escape View File

@ -108,75 +108,75 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
endif
if (prec%iprcparm(smooth_sweeps_) == 1) then
if (prec%iprcparm(mld_smooth_sweeps_) == 1) then
select case(prec%iprcparm(sub_solve_))
case(ilu_n_,milu_n_,ilu_t_)
select case(prec%iprcparm(mld_sub_solve_))
case(mld_ilu_n_,mld_milu_n_,mld_ilu_t_)
select case(toupper(trans))
case('N')
call psb_spsm(zone,prec%av(l_pr_),x,zzero,ww,desc_data,info,&
call psb_spsm(zone,prec%av(mld_l_pr_),x,zzero,ww,desc_data,info,&
& trans='N',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,&
call psb_spsm(alpha,prec%av(mld_u_pr_),ww,beta,y,desc_data,info,&
& trans='N',unit='U',choice=psb_none_, work=aux)
if(info /=0) goto 9999
case('T','C')
call psb_spsm(zone,prec%av(u_pr_),x,zzero,ww,desc_data,info,&
call psb_spsm(zone,prec%av(mld_u_pr_),x,zzero,ww,desc_data,info,&
& 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,&
call psb_spsm(alpha,prec%av(mld_l_pr_),ww,beta,y,desc_data,info,&
& trans=trans,unit='U',choice=psb_none_,work=aux)
if(info /=0) goto 9999
end select
case(slu_)
case(mld_slu_)
ww(1:n_row) = x(1:n_row)
select case(toupper(trans))
case('N')
call mld_zslu_solve(0,n_row,1,ww,n_row,prec%iprcparm(slu_ptr_),info)
call mld_zslu_solve(0,n_row,1,ww,n_row,prec%iprcparm(mld_slu_ptr_),info)
case('T')
call mld_zslu_solve(1,n_row,1,ww,n_row,prec%iprcparm(slu_ptr_),info)
call mld_zslu_solve(1,n_row,1,ww,n_row,prec%iprcparm(mld_slu_ptr_),info)
case('C')
call mld_zslu_solve(2,n_row,1,ww,n_row,prec%iprcparm(slu_ptr_),info)
call mld_zslu_solve(2,n_row,1,ww,n_row,prec%iprcparm(mld_slu_ptr_),info)
end select
if(info /=0) goto 9999
call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
case(sludist_)
case(mld_sludist_)
!!$ write(0,*) 'Calling SLUDist_solve ',n_row
!!$ write(0,*) 'Calling mld_sludist_solve ',n_row
ww(1:n_row) = x(1:n_row)
select case(toupper(trans))
case('N')
call mld_zsludist_solve(0,n_row,1,ww,n_row,prec%iprcparm(slud_ptr_),info)
call mld_zsludist_solve(0,n_row,1,ww,n_row,prec%iprcparm(mld_slud_ptr_),info)
case('T')
call mld_zsludist_solve(1,n_row,1,ww,n_row,prec%iprcparm(slud_ptr_),info)
call mld_zsludist_solve(1,n_row,1,ww,n_row,prec%iprcparm(mld_slud_ptr_),info)
case('C')
call mld_zsludist_solve(2,n_row,1,ww,n_row,prec%iprcparm(slud_ptr_),info)
call mld_zsludist_solve(2,n_row,1,ww,n_row,prec%iprcparm(mld_slud_ptr_),info)
end select
if(info /=0) goto 9999
call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
case (umf_)
case (mld_umf_)
select case(toupper(trans))
case('N')
call mld_zumf_solve(0,n_row,ww,x,n_row,prec%iprcparm(umf_numptr_),info)
call mld_zumf_solve(0,n_row,ww,x,n_row,prec%iprcparm(mld_umf_numptr_),info)
case('T')
call mld_zumf_solve(1,n_row,ww,x,n_row,prec%iprcparm(umf_numptr_),info)
call mld_zumf_solve(1,n_row,ww,x,n_row,prec%iprcparm(mld_umf_numptr_),info)
case('C')
call mld_zumf_solve(2,n_row,ww,x,n_row,prec%iprcparm(umf_numptr_),info)
call mld_zumf_solve(2,n_row,ww,x,n_row,prec%iprcparm(mld_umf_numptr_),info)
end select
if(info /=0) goto 9999
@ -184,15 +184,15 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
case default
write(0,*) 'Unknown factorization type in bjac_aply',prec%iprcparm(sub_solve_)
write(0,*) 'Unknown factorization type in mld_bjac_aply',prec%iprcparm(mld_sub_solve_)
end select
if (debugprt) write(0,*)' Y: ',y(:)
else if (prec%iprcparm(smooth_sweeps_) > 1) then
else if (prec%iprcparm(mld_smooth_sweeps_) > 1) then
! Note: we have to add TRANS to this one !!!!!!!!!
if (size(prec%av) < ap_nd_) then
if (size(prec%av) < mld_ap_nd_) then
info = 4011
goto 9999
endif
@ -207,50 +207,50 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
tx = zzero
ty = zzero
select case(prec%iprcparm(sub_solve_))
case(ilu_n_,milu_n_,ilu_t_)
do i=1, prec%iprcparm(smooth_sweeps_)
select case(prec%iprcparm(mld_sub_solve_))
case(mld_ilu_n_,mld_milu_n_,mld_ilu_t_)
do i=1, prec%iprcparm(mld_smooth_sweeps_)
! X(k+1) = M^-1*(b-N*X(k))
ty(1:n_row) = x(1:n_row)
call psb_spmm(-zone,prec%av(ap_nd_),tx,zone,ty,&
call psb_spmm(-zone,prec%av(mld_ap_nd_),tx,zone,ty,&
& prec%desc_data,info,work=aux)
if(info /=0) goto 9999
call psb_spsm(zone,prec%av(l_pr_),ty,zzero,ww,&
call psb_spsm(zone,prec%av(mld_l_pr_),ty,zzero,ww,&
& prec%desc_data,info,&
& trans='N',unit='L',diag=prec%d,choice=psb_none_,work=aux)
if(info /=0) goto 9999
call psb_spsm(zone,prec%av(u_pr_),ww,zzero,tx,&
call psb_spsm(zone,prec%av(mld_u_pr_),ww,zzero,tx,&
& prec%desc_data,info,&
& trans='N',unit='U',choice=psb_none_,work=aux)
if(info /=0) goto 9999
end do
case(sludist_)
write(0,*) 'No sense in having SLUDist with JAC_SWEEPS >1'
case(mld_sludist_)
write(0,*) 'No sense in having SLUDist with Jmld_ac_SWEEPS >1'
info=4010
goto 9999
case(slu_)
do i=1, prec%iprcparm(smooth_sweeps_)
case(mld_slu_)
do i=1, prec%iprcparm(mld_smooth_sweeps_)
! X(k+1) = M^-1*(b-N*X(k))
ty(1:n_row) = x(1:n_row)
call psb_spmm(-zone,prec%av(ap_nd_),tx,zone,ty,&
call psb_spmm(-zone,prec%av(mld_ap_nd_),tx,zone,ty,&
& prec%desc_data,info,work=aux)
if(info /=0) goto 9999
call mld_zslu_solve(0,n_row,1,ty,n_row,prec%iprcparm(slu_ptr_),info)
call mld_zslu_solve(0,n_row,1,ty,n_row,prec%iprcparm(mld_slu_ptr_),info)
if(info /=0) goto 9999
tx(1:n_row) = ty(1:n_row)
end do
case(umf_)
do i=1, prec%iprcparm(smooth_sweeps_)
case(mld_umf_)
do i=1, prec%iprcparm(mld_smooth_sweeps_)
! X(k+1) = M^-1*(b-N*X(k))
ty(1:n_row) = x(1:n_row)
call psb_spmm(-zone,prec%av(ap_nd_),tx,zone,ty,&
call psb_spmm(-zone,prec%av(mld_ap_nd_),tx,zone,ty,&
& prec%desc_data,info,work=aux)
if(info /=0) goto 9999
call mld_zumf_solve(0,n_row,ww,ty,n_row,&
& prec%iprcparm(umf_numptr_),info)
& prec%iprcparm(mld_umf_numptr_),info)
if(info /=0) goto 9999
tx(1:n_row) = ww(1:n_row)
end do
@ -266,7 +266,7 @@ subroutine mld_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
else
info = 10
call psb_errpush(info,name,&
& i_err=(/2,prec%iprcparm(smooth_sweeps_),0,0,0/))
& i_err=(/2,prec%iprcparm(mld_smooth_sweeps_),0,0,0/))
goto 9999
endif

122
mlprec/mld_zbjac_bld.f90
Unescape Escape View File

@ -94,10 +94,10 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
endif
trans = 'N'
unitd = 'U'
if (p%iprcparm(n_ovr_) < 0) then
if (p%iprcparm(mld_n_ovr_) < 0) then
info = 11
int_err(1) = 1
int_err(2) = p%iprcparm(n_ovr_)
int_err(2) = p%iprcparm(mld_n_ovr_)
call psb_errpush(info,name,i_err=int_err)
goto 9999
endif
@ -108,9 +108,9 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
t1= psb_wtime()
if(debug) write(0,*)me,': calling mld_asmat_bld',&
& p%iprcparm(prec_type_),p%iprcparm(n_ovr_)
& p%iprcparm(mld_prec_type_),p%iprcparm(mld_n_ovr_)
if (debug) call psb_barrier(ictxt)
call mld_asmat_bld(p%iprcparm(prec_type_),p%iprcparm(n_ovr_),a,&
call mld_asmat_bld(p%iprcparm(mld_prec_type_),p%iprcparm(mld_n_ovr_),a,&
& blck,desc_a,upd,p%desc_data,info,outfmt=csrfmt)
if (debugprt) then
@ -132,7 +132,7 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
if (debug) call psb_barrier(ictxt)
select case(p%iprcparm(sub_ren_))
select case(p%iprcparm(mld_sub_ren_))
case (1:)
@ -151,23 +151,23 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
!------------------------------------------------------------------
! Split AC=M+N N off-diagonal part
! Output in COO format.
call psb_sp_clip(atmp,p%av(ap_nd_),info,&
call psb_sp_clip(atmp,p%av(mld_ap_nd_),info,&
& jmin=atmp%m+1,rscale=.false.,cscale=.false.)
call psb_spcnv(p%av(ap_nd_),info,afmt='csr',dupl=psb_dupl_add_)
call psb_spcnv(p%av(mld_ap_nd_),info,afmt='csr',dupl=psb_dupl_add_)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_spcnv csr 1')
goto 9999
end if
k = psb_sp_get_nnzeros(p%av(ap_nd_))
k = psb_sp_get_nnzeros(p%av(mld_ap_nd_))
call psb_sum(ictxt,k)
if (k == 0) then
! If the off diagonal part is emtpy, there's no point
! in doing multiple Jacobi sweeps. This is certain
! to happen when running on a single processor.
p%iprcparm(smooth_sweeps_) = 1
p%iprcparm(mld_smooth_sweeps_) = 1
end if
@ -181,9 +181,9 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
if (debug) write(0,*) me,' Factoring rows ',&
&atmp%m,a%m,blck%m,atmp%ia2(atmp%m+1)-1
select case(p%iprcparm(sub_solve_))
select case(p%iprcparm(mld_sub_solve_))
case(ilu_n_,milu_n_,ilu_t_)
case(mld_ilu_n_,mld_milu_n_,mld_ilu_t_)
call psb_spcnv(atmp,info,afmt='csr',dupl=psb_dupl_add_)
if (info /= 0) then
@ -203,18 +203,18 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
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
call psb_csprt(80+me,p%av(mld_l_pr_),head='% Local L factor')
write(80+me,*) '% Diagonal: ',p%av(mld_l_pr_)%m
do i=1,p%av(mld_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')
call psb_csprt(80+me,p%av(mld_u_pr_),head='% Local U factor')
close(80+me)
endif
case(slu_)
case(mld_slu_)
call psb_spcnv(atmp,info,afmt='csr',dupl=psb_dupl_add_)
if (info /= 0) then
@ -224,11 +224,11 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
call mld_slu_bld(atmp,p%desc_data,p,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='slu_bld')
call psb_errpush(4010,name,a_err='mld_slu_bld')
goto 9999
end if
case(umf_)
case(mld_umf_)
call psb_spcnv(atmp,info,afmt='csc',dupl=psb_dupl_add_)
if (info /= 0) then
@ -237,20 +237,20 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
end if
call mld_umf_bld(atmp,p%desc_data,p,info)
if(debug) write(0,*)me,': Done umf_bld ',info
if(debug) write(0,*)me,': Done mld_umf_bld ',info
if (info /= 0) then
call psb_errpush(4010,name,a_err='umf_bld')
call psb_errpush(4010,name,a_err='mld_umf_bld')
goto 9999
end if
case(f_none_)
case(mld_f_none_)
info=4010
call psb_errpush(info,name,a_err='Inconsistent prec f_none_')
call psb_errpush(info,name,a_err='Inconsistent prec mld_f_none_')
goto 9999
case default
info=4010
call psb_errpush(info,name,a_err='Unknown sub_solve_')
call psb_errpush(info,name,a_err='Unknown mld_sub_solve_')
goto 9999
end select
@ -267,37 +267,37 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
case(0) ! No renumbering
select case(p%iprcparm(sub_solve_))
select case(p%iprcparm(mld_sub_solve_))
case(ilu_n_,milu_n_,ilu_t_)
case(mld_ilu_n_,mld_milu_n_,mld_ilu_t_)
if (p%iprcparm(smooth_sweeps_) > 1) then
if (p%iprcparm(mld_smooth_sweeps_) > 1) then
n_row = psb_cd_get_local_rows(p%desc_data)
n_col = psb_cd_get_local_cols(p%desc_data)
nrow_a = a%m
! The following is known to work
! given that the output from CLIP is in COO.
call psb_sp_clip(a,p%av(ap_nd_),info,&
call psb_sp_clip(a,p%av(mld_ap_nd_),info,&
& jmin=nrow_a+1,rscale=.false.,cscale=.false.)
call psb_sp_clip(blck,atmp,info,&
& jmin=nrow_a+1,rscale=.false.,cscale=.false.)
call psb_rwextd(n_row,p%av(ap_nd_),info,b=atmp)
call psb_rwextd(n_row,p%av(mld_ap_nd_),info,b=atmp)
call psb_spcnv(p%av(ap_nd_),info,afmt='csr',dupl=psb_dupl_add_)
call psb_spcnv(p%av(mld_ap_nd_),info,afmt='csr',dupl=psb_dupl_add_)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_spcnv csr 4')
goto 9999
end if
k = psb_sp_get_nnzeros(p%av(ap_nd_))
k = psb_sp_get_nnzeros(p%av(mld_ap_nd_))
call psb_sum(ictxt,k)
if (k == 0) then
! If the off diagonal part is emtpy, there's no point
! in doing multiple Jacobi sweeps. This is certain
! to happen when running on a single processor.
p%iprcparm(smooth_sweeps_) = 1
p%iprcparm(mld_smooth_sweeps_) = 1
end if
call psb_sp_free(atmp,info)
end if
@ -314,18 +314,18 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
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
call psb_csprt(80+me,p%av(mld_l_pr_),head='% Local L factor')
write(80+me,*) '% Diagonal: ',p%av(mld_l_pr_)%m
do i=1,p%av(mld_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')
call psb_csprt(80+me,p%av(mld_u_pr_),head='% Local U factor')
close(80+me)
endif
case(slu_)
case(mld_slu_)
call psb_spcnv(a,atmp,info,afmt='coo')
if (info /= 0) then
@ -337,34 +337,34 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
n_col = psb_cd_get_local_cols(p%desc_data)
call psb_rwextd(n_row,atmp,info,b=blck)
if (p%iprcparm(smooth_sweeps_) > 1) then
if (p%iprcparm(mld_smooth_sweeps_) > 1) then
!------------------------------------------------------------------
! Split AC=M+N N off-diagonal part
! Output in COO format.
call psb_sp_clip(atmp,p%av(ap_nd_),info,&
call psb_sp_clip(atmp,p%av(mld_ap_nd_),info,&
& jmin=atmp%m+1,rscale=.false.,cscale=.false.)
call psb_spcnv(p%av(ap_nd_),info,afmt='csr',dupl=psb_dupl_add_)
call psb_spcnv(p%av(mld_ap_nd_),info,afmt='csr',dupl=psb_dupl_add_)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_spcnv csr 6')
goto 9999
end if
k = psb_sp_get_nnzeros(p%av(ap_nd_))
k = psb_sp_get_nnzeros(p%av(mld_ap_nd_))
call psb_sum(ictxt,k)
if (k == 0) then
! If the off diagonal part is emtpy, there's no point
! in doing multiple Jacobi sweeps. This is certain
! to happen when running on a single processor.
p%iprcparm(smooth_sweeps_) = 1
p%iprcparm(mld_smooth_sweeps_) = 1
end if
endif
if (info == 0) call psb_spcnv(atmp,info,afmt='csr',dupl=psb_dupl_add_)
if (info == 0) call mld_slu_bld(atmp,p%desc_data,p,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='slu_bld')
call psb_errpush(4010,name,a_err='mld_slu_bld')
goto 9999
end if
@ -375,7 +375,7 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
end if
case(sludist_)
case(mld_sludist_)
call psb_spcnv(a,atmp,info,afmt='coo')
if (info /= 0) then
@ -387,34 +387,34 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
n_col = psb_cd_get_local_cols(p%desc_data)
call psb_rwextd(n_row,atmp,info,b=blck)
if (p%iprcparm(smooth_sweeps_) > 1) then
if (p%iprcparm(mld_smooth_sweeps_) > 1) then
!------------------------------------------------------------------
! Split AC=M+N N off-diagonal part
! Output in COO format.
call psb_sp_clip(atmp,p%av(ap_nd_),info,&
call psb_sp_clip(atmp,p%av(mld_ap_nd_),info,&
& jmin=atmp%m+1,rscale=.false.,cscale=.false.)
call psb_spcnv(p%av(ap_nd_),info,afmt='csr',dupl=psb_dupl_add_)
call psb_spcnv(p%av(mld_ap_nd_),info,afmt='csr',dupl=psb_dupl_add_)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_spcnv csr 7')
goto 9999
end if
k = psb_sp_get_nnzeros(p%av(ap_nd_))
k = psb_sp_get_nnzeros(p%av(mld_ap_nd_))
call psb_sum(ictxt,k)
if (k == 0) then
! If the off diagonal part is emtpy, there's no point
! in doing multiple Jacobi sweeps. This is certain
! to happen when running on a single processor.
p%iprcparm(smooth_sweeps_) = 1
p%iprcparm(mld_smooth_sweeps_) = 1
end if
endif
if (info == 0) call psb_spcnv(atmp,info,afmt='csr',dupl=psb_dupl_add_)
if (info == 0) call mld_sludist_bld(atmp,p%desc_data,p,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='slu_bld')
call psb_errpush(4010,name,a_err='mld_slu_bld')
goto 9999
end if
@ -424,7 +424,7 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
goto 9999
end if
case(umf_)
case(mld_umf_)
call psb_spcnv(a,atmp,info,afmt='coo')
@ -437,28 +437,28 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
n_col = psb_cd_get_local_cols(p%desc_data)
call psb_rwextd(n_row,atmp,info,b=blck)
if (p%iprcparm(smooth_sweeps_) > 1) then
if (p%iprcparm(mld_smooth_sweeps_) > 1) then
!------------------------------------------------------------------
! Split AC=M+N N off-diagonal part
! Output in COO format.
!!$ write(0,*) 'bjac_bld:' size(p%av),ap_nd_
call psb_sp_clip(atmp,p%av(ap_nd_),info,&
!!$ write(0,*) 'mld_bjac_bld:' size(p%av),mld_ap_nd_
call psb_sp_clip(atmp,p%av(mld_ap_nd_),info,&
& jmin=atmp%m+1,rscale=.false.,cscale=.false.)
call psb_spcnv(p%av(ap_nd_),info,afmt='csr',dupl=psb_dupl_add_)
call psb_spcnv(p%av(mld_ap_nd_),info,afmt='csr',dupl=psb_dupl_add_)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_spcnv csr 8')
goto 9999
end if
k = psb_sp_get_nnzeros(p%av(ap_nd_))
k = psb_sp_get_nnzeros(p%av(mld_ap_nd_))
call psb_sum(ictxt,k)
if (k == 0) then
! If the off diagonal part is emtpy, there's no point
! in doing multiple Jacobi sweeps. This is certain
! to happen when running on a single processor.
p%iprcparm(smooth_sweeps_) = 1
p%iprcparm(mld_smooth_sweeps_) = 1
end if
endif
@ -469,9 +469,9 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
end if
call mld_umf_bld(atmp,p%desc_data,p,info)
if(debug) write(0,*)me,': Done umf_bld ',info
if(debug) write(0,*)me,': Done mld_umf_bld ',info
if (info /= 0) then
call psb_errpush(4010,name,a_err='umf_bld')
call psb_errpush(4010,name,a_err='mld_umf_bld')
goto 9999
end if
@ -482,14 +482,14 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info)
end if
case(f_none_)
case(mld_f_none_)
info=4010
call psb_errpush(info,name,a_err='Inconsistent prec f_none_')
call psb_errpush(info,name,a_err='Inconsistent prec mld_f_none_')
goto 9999
case default
info=4010
call psb_errpush(info,name,a_err='Unknown sub_solve_')
call psb_errpush(info,name,a_err='Unknown mld_sub_solve_')
goto 9999
end select

34
mlprec/mld_zilu_bld.f90
Unescape Escape View File

@ -86,7 +86,7 @@ subroutine mld_zilu_bld(a,desc_a,p,upd,info,blck)
unitd = 'U'
if (allocated(p%av)) then
if (size(p%av) < bp_ilu_avsz) then
if (size(p%av) < mld_bp_ilu_avsz_) then
do i=1,size(p%av)
call psb_sp_free(p%av(i),info)
if (info /= 0) then
@ -99,7 +99,7 @@ subroutine mld_zilu_bld(a,desc_a,p,upd,info,blck)
endif
end if
if (.not.allocated(p%av)) then
allocate(p%av(max_avsz),stat=info)
allocate(p%av(mld_max_avsz_),stat=info)
if (info /= 0) then
call psb_errpush(4000,name)
goto 9999
@ -116,12 +116,12 @@ subroutine mld_zilu_bld(a,desc_a,p,upd,info,blck)
if (debug) call psb_barrier(ictxt)
n_row = p%desc_data%matrix_data(psb_n_row_)
p%av(l_pr_)%m = n_row
p%av(l_pr_)%k = n_row
p%av(u_pr_)%m = n_row
p%av(u_pr_)%k = n_row
call psb_sp_all(n_row,n_row,p%av(l_pr_),nztota,info)
if (info == 0) call psb_sp_all(n_row,n_row,p%av(u_pr_),nztota,info)
p%av(mld_l_pr_)%m = n_row
p%av(mld_l_pr_)%k = n_row
p%av(mld_u_pr_)%m = n_row
p%av(mld_u_pr_)%k = n_row
call psb_sp_all(n_row,n_row,p%av(mld_l_pr_),nztota,info)
if (info == 0) call psb_sp_all(n_row,n_row,p%av(mld_u_pr_),nztota,info)
if(info/=0) then
info=4010
ch_err='psb_sp_all'
@ -148,7 +148,7 @@ subroutine mld_zilu_bld(a,desc_a,p,upd,info,blck)
! Ok, factor the matrix.
!
t5 = psb_wtime()
call mld_ilu_fct(p%iprcparm(sub_solve_),a,p%av(l_pr_),p%av(u_pr_),&
call mld_ilu_fct(p%iprcparm(mld_sub_solve_),a,p%av(mld_l_pr_),p%av(mld_u_pr_),&
& p%d,info,blck=blck)
if(info/=0) then
info=4010
@ -164,12 +164,12 @@ subroutine mld_zilu_bld(a,desc_a,p,upd,info,blck)
!
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
call psb_csprt(80+me,p%av(mld_l_pr_),head='% Local L factor')
write(80+me,*) '% Diagonal: ',p%av(mld_l_pr_)%m
do i=1,p%av(mld_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')
call psb_csprt(80+me,p%av(mld_u_pr_),head='% Local U factor')
close(80+me)
endif
@ -183,12 +183,12 @@ subroutine mld_zilu_bld(a,desc_a,p,upd,info,blck)
! 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)
if (psb_sp_getifld(psb_upd_,p%av(mld_u_pr_),info) /= psb_upd_perm_) then
call psb_sp_trim(p%av(mld_u_pr_),info)
endif
if (psb_sp_getifld(psb_upd_,p%av(l_pr_),info) /= psb_upd_perm_) then
call psb_sp_trim(p%av(l_pr_),info)
if (psb_sp_getifld(psb_upd_,p%av(mld_l_pr_),info) /= psb_upd_perm_) then
call psb_sp_trim(p%av(mld_l_pr_),info)
endif

4
mlprec/mld_zilu_fct.f90
Unescape Escape View File

@ -288,7 +288,7 @@ contains
! which means that this entry does not match; thus
! we take it out of diagonal for MILU.
!
if (ialg == milu_n_) then
if (ialg == mld_milu_n_) then
dia = dia - temp*uaspk(jj)
end if
enddo updateloop
@ -428,7 +428,7 @@ contains
! which means that this entry does not match; thus
! we take it out of diagonal for MILU.
!
if (ialg == milu_n_) then
if (ialg == mld_milu_n_) then
dia = dia - temp*uaspk(jj)
end if
enddo updateloopb

116
mlprec/mld_zmlprec_aply.f90
Unescape Escape View File

@ -61,11 +61,11 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
! 1. Number of levels = NLEV = size(baseprecv(:))
! 2. baseprecv(ilev)%av(:) sparse matrices needed for the current level.
! Includes:
! 2.1.: baseprecv(ilev)%av(l_pr_) L factor of ILU preconditioners
! 2.2.: baseprecv(ilev)%av(u_pr_) U factor of ILU preconditioners
! 2.3.: baseprecv(ilev)%av(ap_nd_) Off-diagonal part of A for Jacobi sweeps
! 2.4.: baseprecv(ilev)%av(ac_) Aggregated matrix of level ILEV
! 2.5.: baseprecv(ilev)%av(sm_pr_t_) Smoother prolongator transpose; maps vectors
! 2.1.: baseprecv(ilev)%av(mld_l_pr_) L factor of ILU preconditioners
! 2.2.: baseprecv(ilev)%av(mld_u_pr_) U factor of ILU preconditioners
! 2.3.: baseprecv(ilev)%av(mld_ap_nd_) Off-diagonal part of A for Jacobi sweeps
! 2.4.: baseprecv(ilev)%av(mld_ac_) Aggregated matrix of level ILEV
! 2.5.: baseprecv(ilev)%av(mld_sm_pr_t_) Smoother prolongator transpose; maps vectors
! (ilev-1) ---> (ilev)
! 2.6.: baseprecv(ilev)%av(sm_pr_) Smoother prolongator; maps vectors
! (ilev) ---> (ilev-1)
@ -74,7 +74,7 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
! 3. baseprecv(ilev)%desc_data comm descriptor for level ILEV
! 4. baseprecv(ilev)%base_a Pointer (really a pointer!) to the base matrix
! baseprecv(ilev)%base_desc of the current level, i.e.: if ILEV=1 then A
! else the aggregated matrix av(ac_); so we have
! else the aggregated matrix av(mld_ac_); so we have
! a unified treatment of residuals. Need this to
! avoid passing explicitly matrix A to the
! outer prec. routine
@ -130,15 +130,15 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
end if
select case(baseprecv(2)%iprcparm(ml_type_))
select case(baseprecv(2)%iprcparm(mld_ml_type_))
case(no_ml_)
case(mld_no_ml_)
! Should not really get here.
call psb_errpush(4010,name,a_err='no_ml_ in mlprc_aply?')
call psb_errpush(4010,name,a_err='mld_no_ml_ in mlprc_aply?')
goto 9999
case(add_ml_)
case(mld_add_ml_)
!
@ -191,9 +191,9 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
mlprec_wrk(ilev)%ty(:) = zzero
ismth = baseprecv(ilev)%iprcparm(aggr_kind_)
icm = baseprecv(ilev)%iprcparm(coarse_mat_)
if (ismth /= no_smooth_) then
ismth = baseprecv(ilev)%iprcparm(mld_aggr_kind_)
icm = baseprecv(ilev)%iprcparm(mld_coarse_mat_)
if (ismth /= mld_no_smooth_) then
!
! Smoothed aggregation
!
@ -201,7 +201,7 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
& info,work=work)
if(info /=0) goto 9999
call psb_csmm(zone,baseprecv(ilev)%av(sm_pr_t_),mlprec_wrk(ilev-1)%x2l,&
call psb_csmm(zone,baseprecv(ilev)%av(mld_sm_pr_t_),mlprec_wrk(ilev-1)%x2l,&
& zzero,mlprec_wrk(ilev)%x2l,info)
if(info /=0) goto 9999
@ -217,10 +217,10 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
end if
if (icm ==repl_mat_) Then
if (icm ==mld_repl_mat_) Then
call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nr2l))
else if (icm/= distr_mat_) Then
write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',icm
else if (icm/= mld_distr_mat_) Then
write(0,*) 'Unknown value for baseprecv(2)%iprcparm(mld_coarse_mat_) ',icm
endif
call mld_baseprec_aply(zone,baseprecv(ilev),&
@ -235,10 +235,10 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
n_col = psb_cd_get_local_cols(baseprecv(ilev-1)%desc_data)
nc2l = psb_cd_get_local_cols(baseprecv(ilev)%desc_data)
nr2l = psb_cd_get_local_rows(baseprecv(ilev)%desc_data)
ismth = baseprecv(ilev)%iprcparm(aggr_kind_)
icm = baseprecv(ilev)%iprcparm(coarse_mat_)
ismth = baseprecv(ilev)%iprcparm(mld_aggr_kind_)
icm = baseprecv(ilev)%iprcparm(mld_coarse_mat_)
if (ismth /= no_smooth_) then
if (ismth /= mld_no_smooth_) then
call psb_csmm(zone,baseprecv(ilev)%av(sm_pr_),mlprec_wrk(ilev)%y2l,&
& zone,mlprec_wrk(ilev-1)%y2l,info)
@ -258,16 +258,16 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
if(info /=0) goto 9999
case(mult_ml)
case(mld_mult_ml_)
!
! Multiplicative multilevel
! Pre/post smoothing versions.
!
select case(baseprecv(2)%iprcparm(smooth_pos_))
select case(baseprecv(2)%iprcparm(mld_smooth_pos_))
case(post_smooth_)
case(mld_post_smooth_)
!
! Post smoothing.
@ -310,8 +310,8 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
n_col = psb_cd_get_local_cols(baseprecv(ilev-1)%desc_data)
nc2l = psb_cd_get_local_cols(baseprecv(ilev)%desc_data)
nr2l = psb_cd_get_local_rows(baseprecv(ilev)%desc_data)
ismth = baseprecv(ilev)%iprcparm(aggr_kind_)
icm = baseprecv(ilev)%iprcparm(coarse_mat_)
ismth = baseprecv(ilev)%iprcparm(mld_aggr_kind_)
icm = baseprecv(ilev)%iprcparm(mld_coarse_mat_)
if (debug) write(0,*) me, 'mlpr_aply starting up sweep ',&
& ilev,allocated(baseprecv(ilev)%iprcparm),n_row,n_col,&
@ -330,7 +330,7 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
mlprec_wrk(ilev)%x2l(:) = zzero
mlprec_wrk(ilev)%y2l(:) = zzero
mlprec_wrk(ilev)%tx(:) = zzero
if (ismth /= no_smooth_) then
if (ismth /= mld_no_smooth_) then
!
! Smoothed aggregation
!
@ -340,7 +340,7 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
& baseprecv(ilev-1)%base_desc,info,work=work)
if(info /=0) goto 9999
call psb_csmm(zone,baseprecv(ilev)%av(sm_pr_t_),mlprec_wrk(ilev-1)%x2l, &
call psb_csmm(zone,baseprecv(ilev)%av(mld_sm_pr_t_),mlprec_wrk(ilev-1)%x2l, &
& zzero,mlprec_wrk(ilev)%x2l,info)
if(info /=0) goto 9999
@ -356,18 +356,18 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
end if
if (debug) write(0,*) me, 'mlpr_aply possible sum in up sweep ', &
& ilev,icm,associated(baseprecv(ilev)%base_desc),repl_mat_
& ilev,icm,associated(baseprecv(ilev)%base_desc),mld_repl_mat_
if (debug) write(0,*) me, 'mlpr_aply geaxpby in up sweep X', &
& ilev,associated(baseprecv(ilev)%base_desc),&
& baseprecv(ilev)%base_desc%matrix_data(psb_n_row_),&
& baseprecv(ilev)%base_desc%matrix_data(psb_n_col_),&
& size(mlprec_wrk(ilev)%tx),size(mlprec_wrk(ilev)%x2l)
if (icm == repl_mat_) Then
if (icm == mld_repl_mat_) Then
if (debug) write(0,*) 'Entering psb_sum ',nr2l
call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nr2l))
else if (icm /= distr_mat_) Then
write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ', icm
else if (icm /= mld_distr_mat_) Then
write(0,*) 'Unknown value for baseprecv(2)%iprcparm(mld_coarse_mat_) ', icm
endif
call psb_geaxpby(zone,mlprec_wrk(ilev)%x2l,zzero,mlprec_wrk(ilev)%tx,&
& baseprecv(ilev)%base_desc,info)
@ -383,11 +383,11 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
do ilev=nlev-1, 1, -1
ismth = baseprecv(ilev+1)%iprcparm(aggr_kind_)
ismth = baseprecv(ilev+1)%iprcparm(mld_aggr_kind_)
n_row = psb_cd_get_local_rows(baseprecv(ilev)%base_desc)
if (ismth /= no_smooth_) then
if (ismth == smooth_prol_) &
if (ismth /= mld_no_smooth_) then
if (ismth == mld_smooth_prol_) &
& call psb_halo(mlprec_wrk(ilev+1)%y2l,baseprecv(ilev+1)%desc_data,&
& info,work=work)
call psb_csmm(zone,baseprecv(ilev+1)%av(sm_pr_),mlprec_wrk(ilev+1)%y2l,&
@ -420,7 +420,7 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
if(info /=0) goto 9999
case(pre_smooth_)
case(mld_pre_smooth_)
!
! Pre smoothing.
@ -473,8 +473,8 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
n_col = psb_cd_get_local_cols(baseprecv(ilev-1)%desc_data)
nc2l = psb_cd_get_local_cols(baseprecv(ilev)%desc_data)
nr2l = psb_cd_get_local_rows(baseprecv(ilev)%desc_data)
ismth = baseprecv(ilev)%iprcparm(aggr_kind_)
icm = baseprecv(ilev)%iprcparm(coarse_mat_)
ismth = baseprecv(ilev)%iprcparm(mld_aggr_kind_)
icm = baseprecv(ilev)%iprcparm(mld_coarse_mat_)
allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%y2l(nc2l),&
& mlprec_wrk(ilev)%x2l(nc2l), stat=info)
@ -490,7 +490,7 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
mlprec_wrk(ilev)%tx(:) = zzero
if (ismth /= no_smooth_) then
if (ismth /= mld_no_smooth_) then
!
!Smoothed Aggregation
!
@ -498,7 +498,7 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
& info,work=work)
if(info /=0) goto 9999
call psb_csmm(zone,baseprecv(ilev)%av(sm_pr_t_),mlprec_wrk(ilev-1)%tx,zzero,&
call psb_csmm(zone,baseprecv(ilev)%av(mld_sm_pr_t_),mlprec_wrk(ilev-1)%tx,zzero,&
& mlprec_wrk(ilev)%x2l,info)
if(info /=0) goto 9999
@ -514,10 +514,10 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
end do
end if
if (icm ==repl_mat_) then
if (icm ==mld_repl_mat_) then
call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nr2l))
else if (icm /= distr_mat_) then
write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ', icm
else if (icm /= mld_distr_mat_) then
write(0,*) 'Unknown value for baseprecv(2)%iprcparm(mld_coarse_mat_) ', icm
endif
@ -537,12 +537,12 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
do ilev = nlev-1, 1, -1
ismth = baseprecv(ilev+1)%iprcparm(aggr_kind_)
ismth = baseprecv(ilev+1)%iprcparm(mld_aggr_kind_)
n_row = psb_cd_get_local_rows(baseprecv(ilev)%base_desc)
if (ismth /= no_smooth_) then
if (ismth /= mld_no_smooth_) then
if (ismth == smooth_prol_) &
if (ismth == mld_smooth_prol_) &
& call psb_halo(mlprec_wrk(ilev+1)%y2l,&
& baseprecv(ilev+1)%desc_data,info,work=work)
call psb_csmm(zone,baseprecv(ilev+1)%av(sm_pr_),mlprec_wrk(ilev+1)%y2l,&
@ -568,7 +568,7 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
case(twoside_smooth_)
case(mld_twoside_smooth_)
!
! Symmetrized smoothing.
@ -629,8 +629,8 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
n_col = psb_cd_get_local_cols(baseprecv(ilev-1)%desc_data)
nc2l = psb_cd_get_local_cols(baseprecv(ilev)%desc_data)
nr2l = psb_cd_get_local_rows(baseprecv(ilev)%desc_data)
ismth = baseprecv(ilev)%iprcparm(aggr_kind_)
icm = baseprecv(ilev)%iprcparm(coarse_mat_)
ismth = baseprecv(ilev)%iprcparm(mld_aggr_kind_)
icm = baseprecv(ilev)%iprcparm(mld_coarse_mat_)
allocate(mlprec_wrk(ilev)%ty(nc2l),mlprec_wrk(ilev)%y2l(nc2l),&
& mlprec_wrk(ilev)%x2l(nc2l), stat=info)
@ -647,7 +647,7 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
mlprec_wrk(ilev)%ty(:) = zzero
if (ismth /= no_smooth_) then
if (ismth /= mld_no_smooth_) then
!
!Smoothed Aggregation
!
@ -655,7 +655,7 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
& info,work=work)
if(info /=0) goto 9999
call psb_csmm(zone,baseprecv(ilev)%av(sm_pr_t_),mlprec_wrk(ilev-1)%ty,zzero,&
call psb_csmm(zone,baseprecv(ilev)%av(mld_sm_pr_t_),mlprec_wrk(ilev-1)%ty,zzero,&
& mlprec_wrk(ilev)%x2l,info)
if(info /=0) goto 9999
@ -671,10 +671,10 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
end do
end if
if (icm == repl_mat_) then
if (icm == mld_repl_mat_) then
call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nr2l))
else if (icm /= distr_mat_) then
write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ', icm
else if (icm /= mld_distr_mat_) then
write(0,*) 'Unknown value for baseprecv(2)%iprcparm(mld_coarse_mat_) ', icm
endif
call psb_geaxpby(zone,mlprec_wrk(ilev)%x2l,zzero,mlprec_wrk(ilev)%tx,&
@ -698,11 +698,11 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
do ilev=nlev-1, 1, -1
ismth = baseprecv(ilev+1)%iprcparm(aggr_kind_)
ismth = baseprecv(ilev+1)%iprcparm(mld_aggr_kind_)
n_row = psb_cd_get_local_rows(baseprecv(ilev)%base_desc)
if (ismth /= no_smooth_) then
if (ismth == smooth_prol_) &
if (ismth /= mld_no_smooth_) then
if (ismth == mld_smooth_prol_) &
& call psb_halo(mlprec_wrk(ilev+1)%y2l,baseprecv(ilev+1)%desc_data,&
& info,work=work)
call psb_csmm(zone,baseprecv(ilev+1)%av(sm_pr_),mlprec_wrk(ilev+1)%y2l,&
@ -737,14 +737,14 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
case default
call psb_errpush(4013,name,a_err='wrong smooth_pos',&
& i_Err=(/baseprecv(2)%iprcparm(smooth_pos_),0,0,0,0/))
& i_Err=(/baseprecv(2)%iprcparm(mld_smooth_pos_),0,0,0,0/))
goto 9999
end select
case default
call psb_errpush(4013,name,a_err='wrong mltype',&
& i_Err=(/baseprecv(2)%iprcparm(ml_type_),0,0,0,0/))
& i_Err=(/baseprecv(2)%iprcparm(mld_ml_type_),0,0,0,0/))
goto 9999
end select

42
mlprec/mld_zmlprec_bld.f90
Unescape Escape View File

@ -68,34 +68,34 @@ subroutine mld_zmlprec_bld(a,desc_a,p,info)
call psb_errpush(info,name)
goto 9999
endif
call mld_check_def(p%iprcparm(ml_type_),'Multilevel type',&
& mult_ml,is_legal_ml_type)
call mld_check_def(p%iprcparm(aggr_alg_),'aggregation',&
& dec_aggr_,is_legal_ml_aggr_kind)
call mld_check_def(p%iprcparm(aggr_kind_),'Smoother kind',&
& smooth_prol_,is_legal_ml_smth_kind)
call mld_check_def(p%iprcparm(coarse_mat_),'Coarse matrix',&
& distr_mat_,is_legal_ml_coarse_mat)
call mld_check_def(p%iprcparm(smooth_pos_),'smooth_pos',&
& pre_smooth_,is_legal_ml_smooth_pos)
call mld_check_def(p%iprcparm(mld_ml_type_),'Multilevel type',&
& mld_mult_ml_,is_legal_ml_type)
call mld_check_def(p%iprcparm(mld_aggr_alg_),'aggregation',&
& mld_dec_aggr_,is_legal_ml_aggr_kind)
call mld_check_def(p%iprcparm(mld_aggr_kind_),'Smoother kind',&
& mld_smooth_prol_,is_legal_ml_smth_kind)
call mld_check_def(p%iprcparm(mld_coarse_mat_),'Coarse matrix',&
& mld_distr_mat_,is_legal_ml_coarse_mat)
call mld_check_def(p%iprcparm(mld_smooth_pos_),'smooth_pos',&
& mld_pre_smooth_,is_legal_ml_smooth_pos)
!!$ nullify(p%desc_data)
select case(p%iprcparm(sub_solve_))
case(ilu_n_)
call mld_check_def(p%iprcparm(sub_fill_in_),'Level',0,is_legal_ml_lev)
case(ilu_t_)
call mld_check_def(p%dprcparm(fact_eps_),'Eps',dzero,is_legal_ml_eps)
select case(p%iprcparm(mld_sub_solve_))
case(mld_ilu_n_)
call mld_check_def(p%iprcparm(mld_sub_fill_in_),'Level',0,is_legal_ml_lev)
case(mld_ilu_t_)
call mld_check_def(p%dprcparm(mld_fact_eps_),'Eps',dzero,is_legal_ml_eps)
end select
call mld_check_def(p%dprcparm(aggr_damp_),'omega',dzero,is_legal_omega)
call mld_check_def(p%iprcparm(smooth_sweeps_),'Jacobi sweeps',&
call mld_check_def(p%dprcparm(mld_aggr_damp_),'omega',dzero,is_legal_omega)
call mld_check_def(p%iprcparm(mld_smooth_sweeps_),'Jacobi sweeps',&
& 1,is_legal_jac_sweeps)
! Currently this is ignored by gen_aggrmap, but it could be
! changed in the future. Need to package nlaggr & mlia in a
! private data structure?
call mld_aggrmap_bld(p%iprcparm(aggr_alg_),a,desc_a,p%nlaggr,p%mlia,info)
call mld_aggrmap_bld(p%iprcparm(mld_aggr_alg_),a,desc_a,p%nlaggr,p%mlia,info)
if(info /= 0) then
info=4010
@ -132,11 +132,11 @@ subroutine mld_zmlprec_bld(a,desc_a,p,info)
! We have used a separate ac because:
! 1. We want to reuse the same routines mld_ilu_bld etc.
! 2. We do NOT want to pass an argument twice to them
! p%av(ac_) and p, as this would violate the Fortran standard
! p%av(mld_ac_) and p, as this would violate the Fortran standard
! Hence a separate AC and a TRANSFER function at the end.
!
call psb_sp_transfer(ac,p%av(ac_),info)
p%base_a => p%av(ac_)
call psb_sp_transfer(ac,p%av(mld_ac_),info)
p%base_a => p%av(mld_ac_)
call psb_cdtransfer(desc_ac,p%desc_ac,info)
if (info /= 0) then

2
mlprec/mld_zprecbld.f90
Unescape Escape View File

@ -157,7 +157,7 @@ contains
if (allocated(p%av)) then
! Have not decided what to do yet
end if
allocate(p%av(max_avsz),stat=info)
allocate(p%av(mld_max_avsz_),stat=info)
!!$ if (info /= 0) return
do k=1,size(p%av)
call psb_nullify_sp(p%av(k))

164
mlprec/mld_zprecinit.f90
Unescape Escape View File

@ -61,65 +61,65 @@ subroutine mld_zprecinit(p,ptype,info,nlev)
nlev_ = 1
ilev_ = 1
allocate(p%baseprecv(nlev_),stat=info)
if (info == 0) call psb_realloc(ifpsz,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(dfpsz,p%baseprecv(ilev_)%dprcparm,info)
if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(mld_dfpsz_,p%baseprecv(ilev_)%dprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(prec_type_) = noprec_
p%baseprecv(ilev_)%iprcparm(sub_solve_) = f_none_
p%baseprecv(ilev_)%iprcparm(sub_restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(sub_prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(sub_ren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(smooth_sweeps_) = 1
p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_noprec_
p%baseprecv(ilev_)%iprcparm(mld_sub_solve_) = mld_f_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_ren_) = 0
p%baseprecv(ilev_)%iprcparm(mld_n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(mld_smooth_sweeps_) = 1
case ('DIAG')
nlev_ = 1
ilev_ = 1
allocate(p%baseprecv(nlev_),stat=info)
if (info == 0) call psb_realloc(ifpsz,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(dfpsz,p%baseprecv(ilev_)%dprcparm,info)
if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(mld_dfpsz_,p%baseprecv(ilev_)%dprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(prec_type_) = diag_
p%baseprecv(ilev_)%iprcparm(sub_solve_) = f_none_
p%baseprecv(ilev_)%iprcparm(sub_restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(sub_prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(sub_ren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(smooth_sweeps_) = 1
p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_diag_
p%baseprecv(ilev_)%iprcparm(mld_sub_solve_) = mld_f_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_ren_) = 0
p%baseprecv(ilev_)%iprcparm(mld_n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(mld_smooth_sweeps_) = 1
case ('BJAC')
nlev_ = 1
ilev_ = 1
allocate(p%baseprecv(nlev_),stat=info)
if (info == 0) call psb_realloc(ifpsz,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(dfpsz,p%baseprecv(ilev_)%dprcparm,info)
if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(mld_dfpsz_,p%baseprecv(ilev_)%dprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(prec_type_) = bjac_
p%baseprecv(ilev_)%iprcparm(sub_solve_) = ilu_n_
p%baseprecv(ilev_)%iprcparm(sub_restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(sub_prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(sub_ren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(sub_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(smooth_sweeps_) = 1
p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_bjac_
p%baseprecv(ilev_)%iprcparm(mld_sub_solve_) = mld_ilu_n_
p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_ren_) = 0
p%baseprecv(ilev_)%iprcparm(mld_n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(mld_sub_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(mld_smooth_sweeps_) = 1
case ('ASM','AS')
nlev_ = 1
ilev_ = 1
allocate(p%baseprecv(nlev_),stat=info)
if (info == 0) call psb_realloc(ifpsz,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(dfpsz,p%baseprecv(ilev_)%dprcparm,info)
if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(mld_dfpsz_,p%baseprecv(ilev_)%dprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(prec_type_) = as_
p%baseprecv(ilev_)%iprcparm(sub_solve_) = ilu_n_
p%baseprecv(ilev_)%iprcparm(sub_restr_) = psb_halo_
p%baseprecv(ilev_)%iprcparm(sub_prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(sub_ren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 1
p%baseprecv(ilev_)%iprcparm(sub_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(smooth_sweeps_) = 1
p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_as_
p%baseprecv(ilev_)%iprcparm(mld_sub_solve_) = mld_ilu_n_
p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_halo_
p%baseprecv(ilev_)%iprcparm(mld_sub_prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_ren_) = 0
p%baseprecv(ilev_)%iprcparm(mld_n_ovr_) = 1
p%baseprecv(ilev_)%iprcparm(mld_sub_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(mld_smooth_sweeps_) = 1
case ('MLD', 'ML')
@ -134,58 +134,58 @@ subroutine mld_zprecinit(p,ptype,info,nlev)
endif
ilev_ = 1
allocate(p%baseprecv(nlev_),stat=info)
if (info == 0) call psb_realloc(ifpsz,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(dfpsz,p%baseprecv(ilev_)%dprcparm,info)
if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(mld_dfpsz_,p%baseprecv(ilev_)%dprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(prec_type_) = as_
p%baseprecv(ilev_)%iprcparm(sub_solve_) = ilu_n_
p%baseprecv(ilev_)%iprcparm(sub_restr_) = psb_halo_
p%baseprecv(ilev_)%iprcparm(sub_prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(sub_ren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 1
p%baseprecv(ilev_)%iprcparm(sub_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(smooth_sweeps_) = 1
p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_as_
p%baseprecv(ilev_)%iprcparm(mld_sub_solve_) = mld_ilu_n_
p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_halo_
p%baseprecv(ilev_)%iprcparm(mld_sub_prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_ren_) = 0
p%baseprecv(ilev_)%iprcparm(mld_n_ovr_) = 1
p%baseprecv(ilev_)%iprcparm(mld_sub_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(mld_smooth_sweeps_) = 1
if (nlev_ == 1) return
do ilev_ = 2, nlev_ -1
if (info == 0) call psb_realloc(ifpsz,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(dfpsz,p%baseprecv(ilev_)%dprcparm,info)
if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(mld_dfpsz_,p%baseprecv(ilev_)%dprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(prec_type_) = bjac_
p%baseprecv(ilev_)%iprcparm(sub_restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(sub_prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(sub_ren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(ml_type_) = mult_ml
p%baseprecv(ilev_)%iprcparm(aggr_alg_) = dec_aggr_
p%baseprecv(ilev_)%iprcparm(aggr_kind_) = smooth_prol_
p%baseprecv(ilev_)%iprcparm(coarse_mat_) = distr_mat_
p%baseprecv(ilev_)%iprcparm(smooth_pos_) = post_smooth_
p%baseprecv(ilev_)%iprcparm(aggr_eig_) = max_norm_
p%baseprecv(ilev_)%iprcparm(sub_solve_) = ilu_n_
p%baseprecv(ilev_)%iprcparm(sub_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(smooth_sweeps_) = 1
p%baseprecv(ilev_)%dprcparm(aggr_damp_) = 4.d0/3.d0
p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_bjac_
p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_ren_) = 0
p%baseprecv(ilev_)%iprcparm(mld_n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(mld_ml_type_) = mld_mult_ml_
p%baseprecv(ilev_)%iprcparm(mld_aggr_alg_) = mld_dec_aggr_
p%baseprecv(ilev_)%iprcparm(mld_aggr_kind_) = mld_smooth_prol_
p%baseprecv(ilev_)%iprcparm(mld_coarse_mat_) = mld_distr_mat_
p%baseprecv(ilev_)%iprcparm(mld_smooth_pos_) = mld_post_smooth_
p%baseprecv(ilev_)%iprcparm(mld_aggr_eig_) = mld_max_norm_
p%baseprecv(ilev_)%iprcparm(mld_sub_solve_) = mld_ilu_n_
p%baseprecv(ilev_)%iprcparm(mld_sub_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(mld_smooth_sweeps_) = 1
p%baseprecv(ilev_)%dprcparm(mld_aggr_damp_) = 4.d0/3.d0
end do
ilev_ = nlev_
if (info == 0) call psb_realloc(ifpsz,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(dfpsz,p%baseprecv(ilev_)%dprcparm,info)
if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info)
if (info == 0) call psb_realloc(mld_dfpsz_,p%baseprecv(ilev_)%dprcparm,info)
if (info /= 0) return
p%baseprecv(ilev_)%iprcparm(prec_type_) = bjac_
p%baseprecv(ilev_)%iprcparm(sub_restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(sub_prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(sub_ren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(ml_type_) = mult_ml
p%baseprecv(ilev_)%iprcparm(aggr_alg_) = dec_aggr_
p%baseprecv(ilev_)%iprcparm(aggr_kind_) = smooth_prol_
p%baseprecv(ilev_)%iprcparm(coarse_mat_) = distr_mat_
p%baseprecv(ilev_)%iprcparm(smooth_pos_) = post_smooth_
p%baseprecv(ilev_)%iprcparm(aggr_eig_) = max_norm_
p%baseprecv(ilev_)%iprcparm(sub_solve_) = umf_
p%baseprecv(ilev_)%iprcparm(sub_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(smooth_sweeps_) = 4
p%baseprecv(ilev_)%dprcparm(aggr_damp_) = 4.d0/3.d0
p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_bjac_
p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(mld_sub_ren_) = 0
p%baseprecv(ilev_)%iprcparm(mld_n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(mld_ml_type_) = mld_mult_ml_
p%baseprecv(ilev_)%iprcparm(mld_aggr_alg_) = mld_dec_aggr_
p%baseprecv(ilev_)%iprcparm(mld_aggr_kind_) = mld_smooth_prol_
p%baseprecv(ilev_)%iprcparm(mld_coarse_mat_) = mld_distr_mat_
p%baseprecv(ilev_)%iprcparm(mld_smooth_pos_) = mld_post_smooth_
p%baseprecv(ilev_)%iprcparm(mld_aggr_eig_) = mld_max_norm_
p%baseprecv(ilev_)%iprcparm(mld_sub_solve_) = mld_umf_
p%baseprecv(ilev_)%iprcparm(mld_sub_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(mld_smooth_sweeps_) = 4
p%baseprecv(ilev_)%dprcparm(mld_aggr_damp_) = 4.d0/3.d0
case default
write(0,*) 'Unknown preconditioner type request "',ptype,'"'

42
mlprec/mld_zprecset.f90
Unescape Escape View File

@ -80,7 +80,7 @@ subroutine mld_zprecseti(p,what,val,info,ilev)
if (ilev_ == 1) then
! Rules for fine level are slightly different.
select case(what)
case(prec_type_,sub_solve_,sub_restr_,sub_prol_,sub_ren_,n_ovr_,sub_fill_in_,smooth_sweeps_)
case(mld_prec_type_,mld_sub_solve_,mld_sub_restr_,mld_sub_prol_,mld_sub_ren_,mld_n_ovr_,mld_sub_fill_in_,mld_smooth_sweeps_)
p%baseprecv(ilev_)%iprcparm(what) = val
case default
write(0,*) 'Error: trying to call PRECSET with an invalid WHAT'
@ -88,31 +88,31 @@ subroutine mld_zprecseti(p,what,val,info,ilev)
end select
else if (ilev_ > 1) then
select case(what)
case(prec_type_,sub_solve_,sub_restr_,sub_prol_,sub_ren_,n_ovr_,sub_fill_in_,&
& smooth_sweeps_,ml_type_,aggr_alg_,aggr_kind_,coarse_mat_,&
& smooth_pos_,aggr_eig_)
case(mld_prec_type_,mld_sub_solve_,mld_sub_restr_,mld_sub_prol_,mld_sub_ren_,mld_n_ovr_,mld_sub_fill_in_,&
& mld_smooth_sweeps_,mld_ml_type_,mld_aggr_alg_,mld_aggr_kind_,mld_coarse_mat_,&
& mld_smooth_pos_,mld_aggr_eig_)
p%baseprecv(ilev_)%iprcparm(what) = val
case(coarse_solve_)
case(mld_coarse_solve_)
if (ilev_ /= nlev_) then
write(0,*) 'Inconsistent specification of WHAT vs. ILEV'
info = -2
return
end if
p%baseprecv(ilev_)%iprcparm(sub_solve_) = val
case(coarse_sweeps_)
p%baseprecv(ilev_)%iprcparm(mld_sub_solve_) = val
case(mld_coarse_sweeps_)
if (ilev_ /= nlev_) then
write(0,*) 'Inconsistent specification of WHAT vs. ILEV'
info = -2
return
end if
p%baseprecv(ilev_)%iprcparm(smooth_sweeps_) = val
case(coarse_fill_in_)
p%baseprecv(ilev_)%iprcparm(mld_smooth_sweeps_) = val
case(mld_coarse_fill_in_)
if (ilev_ /= nlev_) then
write(0,*) 'Inconsistent specification of WHAT vs. ILEV'
info = -2
return
end if
p%baseprecv(ilev_)%iprcparm(sub_fill_in_) = val
p%baseprecv(ilev_)%iprcparm(mld_sub_fill_in_) = val
case default
write(0,*) 'Error: trying to call PRECSET with an invalid WHAT'
info = -2
@ -122,9 +122,9 @@ subroutine mld_zprecseti(p,what,val,info,ilev)
else if (.not.present(ilev)) then
select case(what)
case(prec_type_,sub_solve_,sub_restr_,sub_prol_,sub_ren_,n_ovr_,sub_fill_in_,&
& smooth_sweeps_,ml_type_,aggr_alg_,aggr_kind_,coarse_mat_,&
& smooth_pos_,aggr_eig_)
case(mld_prec_type_,mld_sub_solve_,mld_sub_restr_,mld_sub_prol_,mld_sub_ren_,mld_n_ovr_,mld_sub_fill_in_,&
& mld_smooth_sweeps_,mld_ml_type_,mld_aggr_alg_,mld_aggr_kind_,mld_coarse_mat_,&
& mld_smooth_pos_,mld_aggr_eig_)
do ilev_=1,nlev_-1
if (.not.allocated(p%baseprecv(ilev_)%iprcparm)) then
write(0,*) 'Error: trying to call PRECSET on an uninitialized preconditioner'
@ -133,27 +133,27 @@ subroutine mld_zprecseti(p,what,val,info,ilev)
endif
p%baseprecv(ilev_)%iprcparm(what) = val
end do
case(coarse_solve_)
case(mld_coarse_solve_)
if (.not.allocated(p%baseprecv(nlev_)%iprcparm)) then
write(0,*) 'Error: trying to call PRECSET on an uninitialized preconditioner'
info = -1
return
endif
p%baseprecv(nlev_)%iprcparm(sub_solve_) = val
case(coarse_sweeps_)
p%baseprecv(nlev_)%iprcparm(mld_sub_solve_) = val
case(mld_coarse_sweeps_)
if (.not.allocated(p%baseprecv(nlev_)%iprcparm)) then
write(0,*) 'Error: trying to call PRECSET on an uninitialized preconditioner'
info = -1
return
endif
p%baseprecv(nlev_)%iprcparm(smooth_sweeps_) = val
case(coarse_fill_in_)
p%baseprecv(nlev_)%iprcparm(mld_smooth_sweeps_) = val
case(mld_coarse_fill_in_)
if (.not.allocated(p%baseprecv(nlev_)%iprcparm)) then
write(0,*) 'Error: trying to call PRECSET on an uninitialized preconditioner'
info = -1
return
endif
p%baseprecv(nlev_)%iprcparm(sub_fill_in_) = val
p%baseprecv(nlev_)%iprcparm(mld_sub_fill_in_) = val
case default
write(0,*) 'Error: trying to call PRECSET with an invalid WHAT'
info = -2
@ -206,14 +206,14 @@ subroutine mld_zprecsetd(p,what,val,info,ilev)
! Rules for fine level are slightly different.
select case(what)
! Right now we don't have any at base level. Will change when
! we implement ilu_t_
! we implement mld_ilu_t_
case default
write(0,*) 'Error: trying to call PRECSET with an invalid WHAT'
info = -2
end select
else if (ilev_ > 1) then
select case(what)
case(aggr_damp_)
case(mld_aggr_damp_)
p%baseprecv(ilev_)%dprcparm(what) = val
case default
write(0,*) 'Error: trying to call PRECSET with an invalid WHAT'

10
mlprec/mld_zslu_bld.f90
Unescape Escape View File

@ -59,7 +59,7 @@ subroutine mld_zslu_bld(a,desc_a,p,info)
call psb_info(ictxt, me, np)
if (toupper(a%fida) /= 'CSR') then
write(0,*) 'Unimplemented input to SLU_BLD'
write(0,*) 'Unimplemented input to mld_slu_BLD'
goto 9999
endif
@ -67,22 +67,22 @@ subroutine mld_zslu_bld(a,desc_a,p,info)
nzt = psb_sp_get_nnzeros(a)
if (Debug) then
write(0,*) me,'Calling psb_slu_factor ',nzt,a%m,&
write(0,*) me,'Calling psb_mld_slu_factor ',nzt,a%m,&
& a%k,p%desc_data%matrix_data(psb_n_row_)
call psb_barrier(ictxt)
endif
call mld_zslu_factor(a%m,nzt,&
& a%aspk,a%ia2,a%ia1,p%iprcparm(slu_ptr_),info)
& a%aspk,a%ia2,a%ia1,p%iprcparm(mld_slu_ptr_),info)
if (info /= 0) then
ch_err='psb_slu_fact'
ch_err='psb_mld_slu_fact'
call psb_errpush(4110,name,a_err=ch_err,i_err=(/info,0,0,0,0/))
goto 9999
end if
if (Debug) then
write(0,*) me, 'SPLUBLD: Done slu_Factor',info,p%iprcparm(slu_ptr_)
write(0,*) me, 'SPLUBLD: Done mld_slu_Factor',info,p%iprcparm(mld_slu_ptr_)
call psb_barrier(ictxt)
endif

4
mlprec/mld_zslud_bld.f90
Unescape Escape View File

@ -60,7 +60,7 @@ subroutine mld_zsludist_bld(a,desc_a,p,info)
call psb_info(ictxt, me, np)
if (toupper(a%fida) /= 'CSR') then
write(0,*) 'Unimplemented input to SLU_BLD'
write(0,*) 'Unimplemented input to mld_slu_BLD'
goto 9999
endif
@ -97,7 +97,7 @@ subroutine mld_zsludist_bld(a,desc_a,p,info)
end do
!!$ write(0,*) 'Process grid : ',npr,npc
call mld_zsludist_factor(mglob,nrow,nzt,ifrst,&
& a%aspk,a%ia2,a%ia1,p%iprcparm(slud_ptr_),&
& a%aspk,a%ia2,a%ia1,p%iprcparm(mld_slud_ptr_),&
& npr, npc, info)
if (info /= 0) then
ch_err='psb_slud_fact'

6
mlprec/mld_zsp_renum.f90
Unescape Escape View File

@ -76,7 +76,7 @@ subroutine mld_zsp_renum(a,desc_a,blck,p,atmp,info)
call psb_spcnv(a,atmp,info,afmt='coo',dupl=psb_dupl_add_)
call psb_rwextd(a%m+blck%m,atmp,info,blck)
if (p%iprcparm(sub_ren_)==renum_glb_) then
if (p%iprcparm(mld_sub_ren_)==mld_renum_glb_) then
! This is the renumbering coherent with global indices..
mglob = psb_cd_get_global_rows(desc_a)
@ -106,7 +106,7 @@ subroutine mld_zsp_renum(a,desc_a,blck,p,atmp,info)
enddo
t3 = psb_wtime()
else if (p%iprcparm(sub_ren_)==renum_gps_) then
else if (p%iprcparm(mld_sub_ren_)==mld_renum_gps_) then
call psb_spcnv(atmp,info,afmt='csr',dupl=psb_dupl_add_)
nztmp = psb_sp_get_nnzeros(atmp)
@ -140,7 +140,7 @@ subroutine mld_zsp_renum(a,desc_a,blck,p,atmp,info)
itmp(1:8) = 0
! write(0,*) me,' Renumbering: Calling Metis'
! write(0,*) size(p%av(u_pr_)%pl),size(p%av(l_pr_)%pr)
! write(0,*) size(p%av(mld_u_pr_)%pl),size(p%av(mld_l_pr_)%pr)
call gps_reduction(atmp%m,atmp%ia2,atmp%ia1,p%perm,p%invperm,info)
if(info/=0) then
info=4010

10
mlprec/mld_zumf_bld.f90
Unescape Escape View File

@ -59,7 +59,7 @@ subroutine mld_zumf_bld(a,desc_a,p,info)
call psb_info(ictxt, me, np)
if (toupper(a%fida) /= 'CSC') then
write(0,*) 'Unimplemented input to UMF_BLD'
write(0,*) 'Unimplemented input to mld_umf_BLD'
goto 9999
endif
@ -67,7 +67,7 @@ subroutine mld_zumf_bld(a,desc_a,p,info)
nzt = psb_sp_get_nnzeros(a)
if (Debug) then
write(0,*) me,'Calling psb_umf_factor ',nzt,a%m,&
write(0,*) me,'Calling psb_mld_umf_factor ',nzt,a%m,&
& a%k,p%desc_data%matrix_data(psb_n_row_)
open(80+me)
call psb_csprt(80+me,a)
@ -77,17 +77,17 @@ subroutine mld_zumf_bld(a,desc_a,p,info)
call mld_zumf_factor(a%m,nzt,&
& a%aspk,a%ia1,a%ia2,&
& p%iprcparm(umf_symptr_),p%iprcparm(umf_numptr_),info)
& p%iprcparm(mld_umf_symptr_),p%iprcparm(mld_umf_numptr_),info)
if (info /= 0) then
i_err(1) = info
info=4110
call psb_errpush(info,name,a_err='psb_umf_fact',i_err=i_err)
call psb_errpush(info,name,a_err='psb_mld_umf_fact',i_err=i_err)
goto 9999
end if
if (Debug) then
write(0,*) me, 'UMFBLD: Done umf_Factor',info,p%iprcparm(umf_numptr_)
write(0,*) me, 'UMFBLD: Done mld_umf_Factor',info,p%iprcparm(mld_umf_numptr_)
call psb_barrier(ictxt)
endif

24
test/fileread/df_bench.f90
Unescape Escape View File

@ -216,22 +216,22 @@ program df_bench
call mld_precinit(pre,precs(pp)%lv2,info,nlev=nlev)
! Defaults are OK for all intermediate levels. Only fix last level.
if (precs(pp)%omega>=0.0) then
call mld_precset(pre,aggr_damp_,precs(pp)%omega,info,ilev=nlev)
call mld_precset(pre,mld_aggr_damp_,precs(pp)%omega,info,ilev=nlev)
end if
call mld_precset(pre,ml_type_, precs(pp)%mltype, info,ilev=nlev)
call mld_precset(pre,aggr_alg_, precs(pp)%aggr, info,ilev=nlev)
call mld_precset(pre,coarse_mat_, precs(pp)%cmat, info,ilev=nlev)
call mld_precset(pre,smooth_pos_, precs(pp)%smthpos, info,ilev=nlev)
call mld_precset(pre,sub_solve_, precs(pp)%ftype2, info,ilev=nlev)
call mld_precset(pre,smooth_sweeps_, precs(pp)%jswp, info,ilev=nlev)
call mld_precset(pre,aggr_kind_, precs(pp)%smthkind,info,ilev=nlev)
call mld_precset(pre,mld_ml_type_, precs(pp)%mltype, info,ilev=nlev)
call mld_precset(pre,mld_aggr_alg_, precs(pp)%aggr, info,ilev=nlev)
call mld_precset(pre,mld_coarse_mat_, precs(pp)%cmat, info,ilev=nlev)
call mld_precset(pre,mld_smooth_pos_, precs(pp)%smthpos, info,ilev=nlev)
call mld_precset(pre,mld_sub_solve_, precs(pp)%ftype2, info,ilev=nlev)
call mld_precset(pre,mld_smooth_sweeps_, precs(pp)%jswp, info,ilev=nlev)
call mld_precset(pre,mld_aggr_kind_, precs(pp)%smthkind,info,ilev=nlev)
else
call mld_precinit(pre,precs(pp)%lv1,info)
end if
call mld_precset(pre,n_ovr_, precs(pp)%novr,info ,ilev=1)
call mld_precset(pre,sub_restr_, precs(pp)%restr,info ,ilev=1)
call mld_precset(pre,sub_prol_, precs(pp)%prol,info ,ilev=1)
call mld_precset(pre,sub_solve_, precs(pp)%ftype1,info ,ilev=1)
call mld_precset(pre,mld_n_ovr_, precs(pp)%novr,info ,ilev=1)
call mld_precset(pre,mld_sub_restr_, precs(pp)%restr,info ,ilev=1)
call mld_precset(pre,mld_sub_prol_, precs(pp)%prol,info ,ilev=1)
call mld_precset(pre,mld_sub_solve_, precs(pp)%ftype1,info ,ilev=1)
! setting initial guess to zero

24
test/fileread/zf_bench.f90
Unescape Escape View File

@ -201,22 +201,22 @@ program zf_bench
call mld_precinit(pre,precs(pp)%lv2,info,nlev=nlev)
! Defaults are OK for all intermediate levels. Only fix last level.
if (precs(pp)%omega>=0.0) then
call mld_precset(pre,aggr_damp_,precs(pp)%omega,info,ilev=nlev)
call mld_precset(pre,mld_aggr_damp_,precs(pp)%omega,info,ilev=nlev)
end if
call mld_precset(pre,ml_type_, precs(pp)%mltype, info,ilev=nlev)
call mld_precset(pre,aggr_alg_, precs(pp)%aggr, info,ilev=nlev)
call mld_precset(pre,coarse_mat_, precs(pp)%cmat, info,ilev=nlev)
call mld_precset(pre,smooth_pos_, precs(pp)%smthpos, info,ilev=nlev)
call mld_precset(pre,sub_solve_, precs(pp)%ftype2, info,ilev=nlev)
call mld_precset(pre,smooth_sweeps_, precs(pp)%jswp, info,ilev=nlev)
call mld_precset(pre,aggr_kind_, precs(pp)%smthkind,info,ilev=nlev)
call mld_precset(pre,mld_ml_type_, precs(pp)%mltype, info,ilev=nlev)
call mld_precset(pre,mld_aggr_alg_, precs(pp)%aggr, info,ilev=nlev)
call mld_precset(pre,mld_coarse_mat_, precs(pp)%cmat, info,ilev=nlev)
call mld_precset(pre,mld_smooth_pos_, precs(pp)%smthpos, info,ilev=nlev)
call mld_precset(pre,mld_sub_solve_, precs(pp)%ftype2, info,ilev=nlev)
call mld_precset(pre,mld_smooth_sweeps_, precs(pp)%jswp, info,ilev=nlev)
call mld_precset(pre,mld_aggr_kind_, precs(pp)%smthkind,info,ilev=nlev)
else
call mld_precinit(pre,precs(pp)%lv1,info)
end if
call mld_precset(pre,n_ovr_, precs(pp)%novr,info ,ilev=1)
call mld_precset(pre,sub_restr_, precs(pp)%restr,info ,ilev=1)
call mld_precset(pre,sub_prol_, precs(pp)%prol,info ,ilev=1)
call mld_precset(pre,sub_solve_, precs(pp)%ftype1,info ,ilev=1)
call mld_precset(pre,mld_n_ovr_, precs(pp)%novr,info ,ilev=1)
call mld_precset(pre,mld_sub_restr_, precs(pp)%restr,info ,ilev=1)
call mld_precset(pre,mld_sub_prol_, precs(pp)%prol,info ,ilev=1)
call mld_precset(pre,mld_sub_solve_, precs(pp)%ftype1,info ,ilev=1)

Write Preview
Loading…
Cancel
Save