mld2p4-2:

examples/fileread/mld_cexample_1lev.f90
 examples/fileread/mld_cexample_ml.f90
 examples/fileread/mld_dexample_1lev.f90
 examples/fileread/mld_dexample_ml.f90
 examples/fileread/mld_sexample_1lev.f90
 examples/fileread/mld_sexample_ml.f90
 examples/fileread/mld_zexample_1lev.f90
 examples/fileread/mld_zexample_ml.f90
 examples/pdegen/mld_dexample_1lev.f90
 examples/pdegen/mld_dexample_ml.f90
 examples/pdegen/mld_sexample_1lev.f90
 examples/pdegen/mld_sexample_ml.f90

Fixed examples.
stopcriterion
Salvatore Filippone 8 years ago
parent 1371fe19c0
commit f9fa6a2849

@ -37,22 +37,21 @@
! POSSIBILITY OF SUCH DAMAGE. ! POSSIBILITY OF SUCH DAMAGE.
! !
! !
! File: mld_cexample_ml.f90 ! File: mld_cexample_1lev.f90
! !
! This sample program solves a linear system by using BiCGStab preconditioned by ! This sample program solves a linear system by using BiCGStab preconditioned by
! RAS with overlap 2 and ILU(0) on the local blocks, as explained in Section 6.1 ! RAS with overlap 2 and ILU(0) on the local blocks, as explained in Section 5.1
! of the MLD2P4 User's and Reference Guide. ! of the MLD2P4 User's and Reference Guide.
! !
! The matrix and the rhs are read from files (if an rhs is not available, the ! The matrix and the rhs are read from files (if an rhs is not available, the
! unit rhs is set). ! unit rhs is set).
! !
program mld_cexample_ml program mld_cexample_1lev
use psb_base_mod use psb_base_mod
use mld_prec_mod use mld_prec_mod
use psb_krylov_mod use psb_krylov_mod
use psb_util_mod use psb_util_mod
use data_input use data_input
implicit none implicit none
! input parameters ! input parameters
@ -86,10 +85,11 @@ program mld_cexample_ml
integer :: i,info,j,m_problem integer :: i,info,j,m_problem
integer(psb_long_int_k_) :: amatsize, precsize, descsize integer(psb_long_int_k_) :: amatsize, precsize, descsize
integer :: ierr, ircode integer :: ierr, ircode
real(psb_dpk_) :: t1, t2, tprec
real(psb_spk_) :: resmx, resmxp real(psb_spk_) :: resmx, resmxp
real(psb_dpk_) :: t1, t2, tprec
character(len=20) :: name character(len=20) :: name
integer, parameter :: iunit=12 integer, parameter :: iunit=12
type(psb_c_vect_type) :: x_col, r_col
! initialize the parallel environment ! initialize the parallel environment
@ -191,22 +191,28 @@ program mld_cexample_ml
write(*,'(" ")') write(*,'(" ")')
end if end if
! set RAS with overlap 2 and ILU(0) on the local blocks
call mld_precinit(P,'AS',info) ! START SETTING PARAMETER
call mld_precset(P,mld_sub_ovr_,2,info)
! set RAS
call P%init('AS',info)
! set number of overlaps
call P%set('SUB_OVR',2,info)
! build the preconditioner ! build the preconditioner
t1 = psb_wtime() t1 = psb_wtime()
call mld_precbld(A,desc_A,P,info) call P%build(A,desc_A,info)
tprec = psb_wtime()-t1 tprec = psb_wtime()-t1
call psb_amx(ictxt, tprec) call psb_amx(ictxt, tprec)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_precbld') call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_precbld')
goto 9999 goto 9999
end if end if
@ -239,8 +245,7 @@ program mld_cexample_ml
call psb_sum(ictxt,descsize) call psb_sum(ictxt,descsize)
call psb_sum(ictxt,precsize) call psb_sum(ictxt,precsize)
call mld_precdescr(P,info) call P%descr(info)
if (iam == psb_root_) then if (iam == psb_root_) then
write(*,'(" ")') write(*,'(" ")')
write(*,'("Matrix: ",A)')mtrx_file write(*,'("Matrix: ",A)')mtrx_file
@ -258,9 +263,9 @@ program mld_cexample_ml
write(*,'("Total memory occupation for PREC : ",i12)')precsize write(*,'("Total memory occupation for PREC : ",i12)')precsize
end if end if
call psb_gather(x_glob,x,desc_a,info,root=psb_root_) call psb_gather(x_glob,x_col,desc_a,info,root=psb_root_)
if (info == psb_success_) & if (info == psb_success_) &
& call psb_gather(r_glob,r,desc_a,info,root=psb_root_) & call psb_gather(r_glob,r_col,desc_a,info,root=psb_root_)
if (info /= psb_success_) goto 9999 if (info /= psb_success_) goto 9999
if (iam == psb_root_) then if (iam == psb_root_) then
write(0,'(" ")') write(0,'(" ")')
@ -284,7 +289,7 @@ program mld_cexample_ml
call psb_gefree(b, desc_A,info) call psb_gefree(b, desc_A,info)
call psb_gefree(x, desc_A,info) call psb_gefree(x, desc_A,info)
call psb_spfree(A, desc_A,info) call psb_spfree(A, desc_A,info)
call mld_precfree(P,info) call P%free(info)
call psb_cdfree(desc_A,info) call psb_cdfree(desc_A,info)
call psb_exit(ictxt) call psb_exit(ictxt)
stop stop
@ -324,4 +329,4 @@ contains
call psb_bcast(ictxt,tol) call psb_bcast(ictxt,tol)
end subroutine get_parms end subroutine get_parms
end program mld_cexample_ml end program mld_cexample_1lev

@ -42,9 +42,19 @@
! This sample program solves a linear system by using BiCGStab coupled with ! This sample program solves a linear system by using BiCGStab coupled with
! one of the following multi-level preconditioner, as explained in Section 6.1 ! one of the following multi-level preconditioner, as explained in Section 6.1
! of the MLD2P4 User's and Reference Guide: ! of the MLD2P4 User's and Reference Guide:
! - choice = 1, default multi-level Schwarz preconditioner (Sec. 6.1, Fig. 2) !
! - choice = 2, hybrid three-level Schwarz preconditioner (Sec. 6.1, Fig. 3) ! - choice = 1, initialize the default multi-level preconditioner solver, i.e.,
! - choice = 3, additive three-level Schwarz preconditioner (Sec. 6.1, Fig. 4) ! V-cycle with basic smoothed aggregation, 1 hybrid forward/backward
! GS sweep as pre/post-smoother and UMFPACK as coarsest-level
! solver(Sec. 5.1, Fig. 2)
!
! - choice = 2, a V-cycle preconditioner with 1 block-Jacobi sweep
! (with ILU(0) on the blocks) as pre- and post-smoother, and 8 block-Jacobi
! sweeps (with ILU(0) on the blocks) as coarsest-level solver(Sec. 5.1, Fig. 3)
!
! - choice = 3, build a W-cycle preconditioner with 2 Gauss-Seidel sweeps as
! post-smoother (and no pre-smoother), a distributed coarsest
! matrix, and MUMPS as coarsest-level solver (Sec. 5.1, Fig. 4)
! !
! The matrix and the rhs are read from files (if an rhs is not available, the ! The matrix and the rhs are read from files (if an rhs is not available, the
! unit rhs is set). ! unit rhs is set).
@ -88,11 +98,11 @@ program mld_cexample_ml
! other variables ! other variables
integer :: choice integer :: choice
integer :: i,info,j,m_problem integer :: i,info,j,m_problem
integer :: ierr, ircode
integer(psb_long_int_k_) :: amatsize, precsize, descsize integer(psb_long_int_k_) :: amatsize, precsize, descsize
real(psb_dpk_) :: t1, t2, tprec integer :: ierr, ircode
real(psb_spk_) :: resmx, resmxp real(psb_spk_) :: resmx, resmxp
character(len=20) :: name real(psb_dpk_) :: t1, t2, tprec
character(len=20) :: name
integer, parameter :: iunit=12 integer, parameter :: iunit=12
! initialize the parallel environment ! initialize the parallel environment
@ -199,36 +209,37 @@ program mld_cexample_ml
case(1) case(1)
! initialize the default multi-level preconditioner, i.e. hybrid ! initialize the default multi-level preconditioner, i.e. V-cycle
! Schwarz, using RAS (with overlap 1 and ILU(0) on the blocks) ! with basic smoothed aggregation, 1 hybrid forward/backward
! as post-smoother and 4 block-Jacobi sweeps (with UMFPACK LU ! GS sweep as pre/post-smoother and UMFPACK as coarsest-level
! on the blocks) as distributed coarse-level solver ! solver
call mld_precinit(P,'ML',info) call P%init('ML',info)
case(2) case(2)
! set a three-level hybrid Schwarz preconditioner, which uses ! initialize a V-cycle preconditioner with 1 block-Jacobi sweep (with
! block Jacobi (with ILU(0) on the blocks) as post-smoother, ! ILU(0) on the blocks) as pre- and post-smoother, and 8 block-Jacobi
! a coarsest matrix replicated on the processors, and the ! sweeps (with ILU(0) on the blocks) as coarsest-level solver
! LU factorization from UMFPACK as coarse-level solver
call mld_precinit(P,'ML',info,nlev=3) call P%init('ML',info)
call mld_precset(P,mld_smoother_type_,'BJAC',info) call P%set('SMOOTHER_TYPE','BJAC',info)
call mld_precset(P,mld_coarse_mat_,'REPL',info) call P%set('COARSE_SOLVE','BJAC',info)
call mld_precset(P,mld_coarse_solve_,'UMF',info) call P%set('COARSE_SWEEPS',8,info)
case(3) case(3)
! set a three-level additive Schwarz preconditioner, which uses ! initialize a W-cycle preconditioner with 2 Gauss-Seidel sweeps as
! RAS (with overlap 1 and ILU(0) on the blocks) as pre- and ! post-smoother (and no pre-smoother), a distributed coarsest
! post-smoother, and 5 block-Jacobi sweeps (with UMFPACK LU ! matrix, and MUMPS as coarsest-level solver
! on the blocks) as distributed coarsest-level solver
call mld_precinit(P,'ML',info,nlev=3) call P%init('ML',info)
call mld_precset(P,mld_ml_type_,'ADD',info) call P%set('ML_TYPE','WCYCLE',info)
call mld_precset(P,mld_smoother_pos_,'TWOSIDE',info) call P%set('SMOOTHER_TYPE','GS',info)
call mld_precset(P,mld_coarse_sweeps_,5,info) call P%set('SMOOTHER_SWEEPS',0,info,pos='PRE')
call P%set('SMOOTHER_SWEEPS',2,info,pos='POST')
call P%set('COARSE_SOLVE','MUMPS',info)
call P%set('COARSE_MAT','DIST',info)
end select end select
@ -237,18 +248,19 @@ program mld_cexample_ml
call psb_barrier(ictxt) call psb_barrier(ictxt)
t1 = psb_wtime() t1 = psb_wtime()
call mld_precbld(A,desc_A,P,info) ! build the preconditioner
call P%hierarchy_build(A,desc_A,info)
call P%smoothers_build(A,desc_A,info)
tprec = psb_wtime()-t1 tprec = psb_wtime()-t1
call psb_amx(ictxt, tprec) call psb_amx(ictxt, tprec)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_precbld') call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_precbld')
goto 9999 goto 9999
end if end if
! set the initial guess ! set the initial guess
call psb_geall(x,desc_A,info) call psb_geall(x,desc_A,info)
call x%zero() call x%zero()
call psb_geasb(x,desc_A,info) call psb_geasb(x,desc_A,info)
@ -276,7 +288,7 @@ program mld_cexample_ml
call psb_sum(ictxt,descsize) call psb_sum(ictxt,descsize)
call psb_sum(ictxt,precsize) call psb_sum(ictxt,precsize)
call mld_precdescr(P,info) call P%descr(info)
if (iam == psb_root_) then if (iam == psb_root_) then
write(*,'(" ")') write(*,'(" ")')
@ -321,7 +333,7 @@ program mld_cexample_ml
call psb_gefree(b, desc_A,info) call psb_gefree(b, desc_A,info)
call psb_gefree(x, desc_A,info) call psb_gefree(x, desc_A,info)
call psb_spfree(A, desc_A,info) call psb_spfree(A, desc_A,info)
call mld_precfree(P,info) call P%free(info)
call psb_cdfree(desc_A,info) call psb_cdfree(desc_A,info)
call psb_exit(ictxt) call psb_exit(ictxt)
stop stop

@ -37,22 +37,21 @@
! POSSIBILITY OF SUCH DAMAGE. ! POSSIBILITY OF SUCH DAMAGE.
! !
! !
! File: mld_dexample_ml.f90 ! File: mld_dexample_1lev.f90
! !
! This sample program solves a linear system by using BiCGStab preconditioned by ! This sample program solves a linear system by using BiCGStab preconditioned by
! RAS with overlap 2 and ILU(0) on the local blocks, as explained in Section 6.1 ! RAS with overlap 2 and ILU(0) on the local blocks, as explained in Section 5.1
! of the MLD2P4 User's and Reference Guide. ! of the MLD2P4 User's and Reference Guide.
! !
! The matrix and the rhs are read from files (if an rhs is not available, the ! The matrix and the rhs are read from files (if an rhs is not available, the
! unit rhs is set). ! unit rhs is set).
! !
program mld_dexample_ml program mld_dexample_1lev
use psb_base_mod use psb_base_mod
use mld_prec_mod use mld_prec_mod
use psb_krylov_mod use psb_krylov_mod
use psb_util_mod use psb_util_mod
use data_input use data_input
use mld_d_mumps_solver
implicit none implicit none
! input parameters ! input parameters
@ -78,7 +77,6 @@ program mld_dexample_ml
real(psb_dpk_) :: tol, err real(psb_dpk_) :: tol, err
integer :: itmax, iter, istop integer :: itmax, iter, istop
integer :: nlev integer :: nlev
type(mld_d_mumps_solver_type) :: mumps_sv
! parallel environment parameters ! parallel environment parameters
integer :: ictxt, iam, np integer :: ictxt, iam, np
@ -87,7 +85,8 @@ program mld_dexample_ml
integer :: i,info,j,m_problem integer :: i,info,j,m_problem
integer(psb_long_int_k_) :: amatsize, precsize, descsize integer(psb_long_int_k_) :: amatsize, precsize, descsize
integer :: ierr, ircode integer :: ierr, ircode
real(psb_dpk_) :: t1, t2, tprec, resmx, resmxp real(psb_dpk_) :: resmx, resmxp
real(psb_dpk_) :: t1, t2, tprec
character(len=20) :: name character(len=20) :: name
integer, parameter :: iunit=12 integer, parameter :: iunit=12
type(psb_d_vect_type) :: x_col, r_col type(psb_d_vect_type) :: x_col, r_col
@ -195,25 +194,25 @@ program mld_dexample_ml
! START SETTING PARAMETER ! START SETTING PARAMETER
! set JAC ! set RAS
call mld_precinit(P,'JAC',info) call P%init('AS',info)
! set MUMPS ad solver
call P%set(mumps_sv,info) ! set number of overlaps
call P%set('SUB_OVR',2,info)
! build the preconditioner ! build the preconditioner
t1 = psb_wtime() t1 = psb_wtime()
call mld_precbld(A,desc_A,P,info) call P%build(A,desc_A,info)
tprec = psb_wtime()-t1 tprec = psb_wtime()-t1
call psb_amx(ictxt, tprec) call psb_amx(ictxt, tprec)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_precbld') call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_precbld')
goto 9999 goto 9999
end if end if
@ -246,8 +245,7 @@ program mld_dexample_ml
call psb_sum(ictxt,descsize) call psb_sum(ictxt,descsize)
call psb_sum(ictxt,precsize) call psb_sum(ictxt,precsize)
call mld_precdescr(P,info) call P%descr(info)
if (iam == psb_root_) then if (iam == psb_root_) then
write(*,'(" ")') write(*,'(" ")')
write(*,'("Matrix: ",A)')mtrx_file write(*,'("Matrix: ",A)')mtrx_file
@ -291,7 +289,7 @@ program mld_dexample_ml
call psb_gefree(b, desc_A,info) call psb_gefree(b, desc_A,info)
call psb_gefree(x, desc_A,info) call psb_gefree(x, desc_A,info)
call psb_spfree(A, desc_A,info) call psb_spfree(A, desc_A,info)
call mld_precfree(P,info) call P%free(info)
call psb_cdfree(desc_A,info) call psb_cdfree(desc_A,info)
call psb_exit(ictxt) call psb_exit(ictxt)
stop stop
@ -331,4 +329,4 @@ contains
call psb_bcast(ictxt,tol) call psb_bcast(ictxt,tol)
end subroutine get_parms end subroutine get_parms
end program mld_dexample_ml end program mld_dexample_1lev

@ -42,9 +42,19 @@
! This sample program solves a linear system by using BiCGStab coupled with ! This sample program solves a linear system by using BiCGStab coupled with
! one of the following multi-level preconditioner, as explained in Section 6.1 ! one of the following multi-level preconditioner, as explained in Section 6.1
! of the MLD2P4 User's and Reference Guide: ! of the MLD2P4 User's and Reference Guide:
! - choice = 1, default multi-level Schwarz preconditioner (Sec. 6.1, Fig. 2) !
! - choice = 2, hybrid three-level Schwarz preconditioner (Sec. 6.1, Fig. 3) ! - choice = 1, initialize the default multi-level preconditioner solver, i.e.,
! - choice = 3, additive three-level Schwarz preconditioner (Sec. 6.1, Fig. 4) ! V-cycle with basic smoothed aggregation, 1 hybrid forward/backward
! GS sweep as pre/post-smoother and UMFPACK as coarsest-level
! solver(Sec. 5.1, Fig. 2)
!
! - choice = 2, a V-cycle preconditioner with 1 block-Jacobi sweep
! (with ILU(0) on the blocks) as pre- and post-smoother, and 8 block-Jacobi
! sweeps (with ILU(0) on the blocks) as coarsest-level solver(Sec. 5.1, Fig. 3)
!
! - choice = 3, build a W-cycle preconditioner with 2 Gauss-Seidel sweeps as
! post-smoother (and no pre-smoother), a distributed coarsest
! matrix, and MUMPS as coarsest-level solver (Sec. 5.1, Fig. 4)
! !
! The matrix and the rhs are read from files (if an rhs is not available, the ! The matrix and the rhs are read from files (if an rhs is not available, the
! unit rhs is set). ! unit rhs is set).
@ -90,7 +100,8 @@ program mld_dexample_ml
integer :: i,info,j,m_problem integer :: i,info,j,m_problem
integer(psb_long_int_k_) :: amatsize, precsize, descsize integer(psb_long_int_k_) :: amatsize, precsize, descsize
integer :: ierr, ircode integer :: ierr, ircode
real(psb_dpk_) :: t1, t2, tprec, resmx, resmxp real(psb_dpk_) :: resmx, resmxp
real(psb_dpk_) :: t1, t2, tprec
character(len=20) :: name character(len=20) :: name
integer, parameter :: iunit=12 integer, parameter :: iunit=12
@ -198,48 +209,37 @@ program mld_dexample_ml
case(1) case(1)
! initialize the default multi-level preconditioner, i.e. hybrid ! initialize the default multi-level preconditioner, i.e. V-cycle
! Schwarz, using RAS (with overlap 1 and ILU(0) on the blocks) ! with basic smoothed aggregation, 1 hybrid forward/backward
! as post-smoother and 4 block-Jacobi sweeps (with UMFPACK LU ! GS sweep as pre/post-smoother and UMFPACK as coarsest-level
! on the blocks) as distributed coarse-level solver ! solver
call mld_precinit(P,'ML',info) call P%init('ML',info)
case(2) case(2)
! set a three-level hybrid Schwarz preconditioner, which uses ! initialize a V-cycle preconditioner with 1 block-Jacobi sweep (with
! block Jacobi (with ILU(0) on the blocks) as post-smoother, ! ILU(0) on the blocks) as pre- and post-smoother, and 8 block-Jacobi
! a coarsest matrix replicated on the processors, and the ! sweeps (with ILU(0) on the blocks) as coarsest-level solver
! LU factorization from UMFPACK as coarse-level solver
call mld_precinit(P,'ML',info,nlev=3) call P%init('ML',info)
call mld_precset(P,mld_smoother_type_,'BJAC',info) call P%set('SMOOTHER_TYPE','BJAC',info)
call mld_precset(P,mld_coarse_mat_,'REPL',info) call P%set('COARSE_SOLVE','BJAC',info)
call mld_precset(P,mld_coarse_solve_,'UMF',info) call P%set('COARSE_SWEEPS',8,info)
case(3) case(3)
! set a three-level additive Schwarz preconditioner, which uses ! initialize a W-cycle preconditioner with 2 Gauss-Seidel sweeps as
! RAS (with overlap 1 and ILU(0) on the blocks) as pre- and ! post-smoother (and no pre-smoother), a distributed coarsest
! post-smoother, and 5 block-Jacobi sweeps (with UMFPACK LU ! matrix, and MUMPS as coarsest-level solver
! on the blocks) as distributed coarsest-level solver
call mld_precinit(P,'ML',info,nlev=3)
call mld_precset(P,mld_ml_type_,'ADD',info)
call mld_precset(P,mld_smoother_pos_,'TWOSIDE',info)
call mld_precset(P,mld_coarse_sweeps_,5,info)
case(4) call P%init('ML',info)
call P%set('ML_TYPE','WCYCLE',info)
! set a three-level hybrid Schwarz preconditioner, which uses call P%set('SMOOTHER_TYPE','GS',info)
! block Jacobi (with ILU(0) on the blocks) as post-smoother, call P%set('SMOOTHER_SWEEPS',0,info,pos='PRE')
! a coarsest matrix replicated on the processors, and the call P%set('SMOOTHER_SWEEPS',2,info,pos='POST')
! multifrontal solver in MUMPS as coarse-level solver call P%set('COARSE_SOLVE','MUMPS',info)
call P%set('COARSE_MAT','DIST',info)
call mld_precinit(P,'ML',info,nlev=3)
call mld_precset(P,mld_smoother_type_,'BJAC',info)
call mld_precset(P,mld_coarse_mat_,'REPL',info)
call mld_precset(P,mld_coarse_solve_,'MUMPS',info)
end select end select
@ -248,13 +248,15 @@ program mld_dexample_ml
call psb_barrier(ictxt) call psb_barrier(ictxt)
t1 = psb_wtime() t1 = psb_wtime()
call mld_precbld(A,desc_A,P,info) ! build the preconditioner
call P%hierarchy_build(A,desc_A,info)
call P%smoothers_build(A,desc_A,info)
tprec = psb_wtime()-t1 tprec = psb_wtime()-t1
call psb_amx(ictxt, tprec) call psb_amx(ictxt, tprec)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_precbld') call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_precbld')
goto 9999 goto 9999
end if end if
@ -286,7 +288,7 @@ program mld_dexample_ml
call psb_sum(ictxt,descsize) call psb_sum(ictxt,descsize)
call psb_sum(ictxt,precsize) call psb_sum(ictxt,precsize)
call mld_precdescr(P,info) call P%descr(info)
if (iam == psb_root_) then if (iam == psb_root_) then
write(*,'(" ")') write(*,'(" ")')
@ -331,7 +333,7 @@ program mld_dexample_ml
call psb_gefree(b, desc_A,info) call psb_gefree(b, desc_A,info)
call psb_gefree(x, desc_A,info) call psb_gefree(x, desc_A,info)
call psb_spfree(A, desc_A,info) call psb_spfree(A, desc_A,info)
call mld_precfree(P,info) call P%free(info)
call psb_cdfree(desc_A,info) call psb_cdfree(desc_A,info)
call psb_exit(ictxt) call psb_exit(ictxt)
stop stop

@ -37,22 +37,21 @@
! POSSIBILITY OF SUCH DAMAGE. ! POSSIBILITY OF SUCH DAMAGE.
! !
! !
! File: mld_sexample_ml.f90 ! File: mld_sexample_1lev.f90
! !
! This sample program solves a linear system by using BiCGStab preconditioned by ! This sample program solves a linear system by using BiCGStab preconditioned by
! RAS with overlap 2 and ILU(0) on the local blocks, as explained in Section 6.1 ! RAS with overlap 2 and ILU(0) on the local blocks, as explained in Section 5.1
! of the MLD2P4 User's and Reference Guide. ! of the MLD2P4 User's and Reference Guide.
! !
! The matrix and the rhs are read from files (if an rhs is not available, the ! The matrix and the rhs are read from files (if an rhs is not available, the
! unit rhs is set). ! unit rhs is set).
! !
program mld_sexample_ml program mld_sexample_1lev
use psb_base_mod use psb_base_mod
use mld_prec_mod use mld_prec_mod
use psb_krylov_mod use psb_krylov_mod
use psb_util_mod use psb_util_mod
use data_input use data_input
implicit none implicit none
! input parameters ! input parameters
@ -86,10 +85,11 @@ program mld_sexample_ml
integer :: i,info,j,m_problem integer :: i,info,j,m_problem
integer(psb_long_int_k_) :: amatsize, precsize, descsize integer(psb_long_int_k_) :: amatsize, precsize, descsize
integer :: ierr, ircode integer :: ierr, ircode
real(psb_dpk_) :: t1, t2, tprec
real(psb_spk_) :: resmx, resmxp real(psb_spk_) :: resmx, resmxp
real(psb_dpk_) :: t1, t2, tprec
character(len=20) :: name character(len=20) :: name
integer, parameter :: iunit=12 integer, parameter :: iunit=12
type(psb_s_vect_type) :: x_col, r_col
! initialize the parallel environment ! initialize the parallel environment
@ -191,22 +191,28 @@ program mld_sexample_ml
write(*,'(" ")') write(*,'(" ")')
end if end if
! set RAS with overlap 2 and ILU(0) on the local blocks
call mld_precinit(P,'AS',info) ! START SETTING PARAMETER
call mld_precset(P,mld_sub_ovr_,2,info)
! set RAS
call P%init('AS',info)
! set number of overlaps
call P%set('SUB_OVR',2,info)
! build the preconditioner ! build the preconditioner
t1 = psb_wtime() t1 = psb_wtime()
call mld_precbld(A,desc_A,P,info) call P%build(A,desc_A,info)
tprec = psb_wtime()-t1 tprec = psb_wtime()-t1
call psb_amx(ictxt, tprec) call psb_amx(ictxt, tprec)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_precbld') call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_precbld')
goto 9999 goto 9999
end if end if
@ -239,8 +245,7 @@ program mld_sexample_ml
call psb_sum(ictxt,descsize) call psb_sum(ictxt,descsize)
call psb_sum(ictxt,precsize) call psb_sum(ictxt,precsize)
call mld_precdescr(P,info) call P%descr(info)
if (iam == psb_root_) then if (iam == psb_root_) then
write(*,'(" ")') write(*,'(" ")')
write(*,'("Matrix: ",A)')mtrx_file write(*,'("Matrix: ",A)')mtrx_file
@ -258,9 +263,9 @@ program mld_sexample_ml
write(*,'("Total memory occupation for PREC : ",i12)')precsize write(*,'("Total memory occupation for PREC : ",i12)')precsize
end if end if
call psb_gather(x_glob,x,desc_a,info,root=psb_root_) call psb_gather(x_glob,x_col,desc_a,info,root=psb_root_)
if (info == psb_success_) & if (info == psb_success_) &
& call psb_gather(r_glob,r,desc_a,info,root=psb_root_) & call psb_gather(r_glob,r_col,desc_a,info,root=psb_root_)
if (info /= psb_success_) goto 9999 if (info /= psb_success_) goto 9999
if (iam == psb_root_) then if (iam == psb_root_) then
write(0,'(" ")') write(0,'(" ")')
@ -284,7 +289,7 @@ program mld_sexample_ml
call psb_gefree(b, desc_A,info) call psb_gefree(b, desc_A,info)
call psb_gefree(x, desc_A,info) call psb_gefree(x, desc_A,info)
call psb_spfree(A, desc_A,info) call psb_spfree(A, desc_A,info)
call mld_precfree(P,info) call P%free(info)
call psb_cdfree(desc_A,info) call psb_cdfree(desc_A,info)
call psb_exit(ictxt) call psb_exit(ictxt)
stop stop
@ -324,4 +329,4 @@ contains
call psb_bcast(ictxt,tol) call psb_bcast(ictxt,tol)
end subroutine get_parms end subroutine get_parms
end program mld_sexample_ml end program mld_sexample_1lev

@ -42,9 +42,19 @@
! This sample program solves a linear system by using BiCGStab coupled with ! This sample program solves a linear system by using BiCGStab coupled with
! one of the following multi-level preconditioner, as explained in Section 6.1 ! one of the following multi-level preconditioner, as explained in Section 6.1
! of the MLD2P4 User's and Reference Guide: ! of the MLD2P4 User's and Reference Guide:
! - choice = 1, default multi-level Schwarz preconditioner (Sec. 6.1, Fig. 2) !
! - choice = 2, hybrid three-level Schwarz preconditioner (Sec. 6.1, Fig. 3) ! - choice = 1, initialize the default multi-level preconditioner solver, i.e.,
! - choice = 3, additive three-level Schwarz preconditioner (Sec. 6.1, Fig. 4) ! V-cycle with basic smoothed aggregation, 1 hybrid forward/backward
! GS sweep as pre/post-smoother and UMFPACK as coarsest-level
! solver(Sec. 5.1, Fig. 2)
!
! - choice = 2, a V-cycle preconditioner with 1 block-Jacobi sweep
! (with ILU(0) on the blocks) as pre- and post-smoother, and 8 block-Jacobi
! sweeps (with ILU(0) on the blocks) as coarsest-level solver(Sec. 5.1, Fig. 3)
!
! - choice = 3, build a W-cycle preconditioner with 2 Gauss-Seidel sweeps as
! post-smoother (and no pre-smoother), a distributed coarsest
! matrix, and MUMPS as coarsest-level solver (Sec. 5.1, Fig. 4)
! !
! The matrix and the rhs are read from files (if an rhs is not available, the ! The matrix and the rhs are read from files (if an rhs is not available, the
! unit rhs is set). ! unit rhs is set).
@ -90,8 +100,8 @@ program mld_sexample_ml
integer :: i,info,j,m_problem integer :: i,info,j,m_problem
integer(psb_long_int_k_) :: amatsize, precsize, descsize integer(psb_long_int_k_) :: amatsize, precsize, descsize
integer :: ierr, ircode integer :: ierr, ircode
real(psb_dpk_) :: t1, t2, tprec real(psb_spk_) :: resmx, resmxp
real(psb_spk_) :: resmx, resmxp real(psb_dpk_) :: t1, t2, tprec
character(len=20) :: name character(len=20) :: name
integer, parameter :: iunit=12 integer, parameter :: iunit=12
@ -199,36 +209,37 @@ program mld_sexample_ml
case(1) case(1)
! initialize the default multi-level preconditioner, i.e. hybrid ! initialize the default multi-level preconditioner, i.e. V-cycle
! Schwarz, using RAS (with overlap 1 and ILU(0) on the blocks) ! with basic smoothed aggregation, 1 hybrid forward/backward
! as post-smoother and 4 block-Jacobi sweeps (with UMFPACK LU ! GS sweep as pre/post-smoother and UMFPACK as coarsest-level
! on the blocks) as distributed coarse-level solver ! solver
call mld_precinit(P,'ML',info) call P%init('ML',info)
case(2) case(2)
! set a three-level hybrid Schwarz preconditioner, which uses ! initialize a V-cycle preconditioner with 1 block-Jacobi sweep (with
! block Jacobi (with ILU(0) on the blocks) as post-smoother, ! ILU(0) on the blocks) as pre- and post-smoother, and 8 block-Jacobi
! a coarsest matrix replicated on the processors, and the ! sweeps (with ILU(0) on the blocks) as coarsest-level solver
! LU factorization from UMFPACK as coarse-level solver
call mld_precinit(P,'ML',info,nlev=3) call P%init('ML',info)
call mld_precset(P,mld_smoother_type_,'BJAC',info) call P%set('SMOOTHER_TYPE','BJAC',info)
call mld_precset(P,mld_coarse_mat_,'REPL',info) call P%set('COARSE_SOLVE','BJAC',info)
call mld_precset(P,mld_coarse_solve_,'UMF',info) call P%set('COARSE_SWEEPS',8,info)
case(3) case(3)
! set a three-level additive Schwarz preconditioner, which uses ! initialize a W-cycle preconditioner with 2 Gauss-Seidel sweeps as
! RAS (with overlap 1 and ILU(0) on the blocks) as pre- and ! post-smoother (and no pre-smoother), a distributed coarsest
! post-smoother, and 5 block-Jacobi sweeps (with UMFPACK LU ! matrix, and MUMPS as coarsest-level solver
! on the blocks) as distributed coarsest-level solver
call mld_precinit(P,'ML',info,nlev=3) call P%init('ML',info)
call mld_precset(P,mld_ml_type_,'ADD',info) call P%set('ML_TYPE','WCYCLE',info)
call mld_precset(P,mld_smoother_pos_,'TWOSIDE',info) call P%set('SMOOTHER_TYPE','GS',info)
call mld_precset(P,mld_coarse_sweeps_,5,info) call P%set('SMOOTHER_SWEEPS',0,info,pos='PRE')
call P%set('SMOOTHER_SWEEPS',2,info,pos='POST')
call P%set('COARSE_SOLVE','MUMPS',info)
call P%set('COARSE_MAT','DIST',info)
end select end select
@ -237,18 +248,19 @@ program mld_sexample_ml
call psb_barrier(ictxt) call psb_barrier(ictxt)
t1 = psb_wtime() t1 = psb_wtime()
call mld_precbld(A,desc_A,P,info) ! build the preconditioner
call P%hierarchy_build(A,desc_A,info)
call P%smoothers_build(A,desc_A,info)
tprec = psb_wtime()-t1 tprec = psb_wtime()-t1
call psb_amx(ictxt, tprec) call psb_amx(ictxt, tprec)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_precbld') call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_precbld')
goto 9999 goto 9999
end if end if
! set the initial guess ! set the initial guess
call psb_geall(x,desc_A,info) call psb_geall(x,desc_A,info)
call x%zero() call x%zero()
call psb_geasb(x,desc_A,info) call psb_geasb(x,desc_A,info)
@ -276,7 +288,7 @@ program mld_sexample_ml
call psb_sum(ictxt,descsize) call psb_sum(ictxt,descsize)
call psb_sum(ictxt,precsize) call psb_sum(ictxt,precsize)
call mld_precdescr(P,info) call P%descr(info)
if (iam == psb_root_) then if (iam == psb_root_) then
write(*,'(" ")') write(*,'(" ")')
@ -321,7 +333,7 @@ program mld_sexample_ml
call psb_gefree(b, desc_A,info) call psb_gefree(b, desc_A,info)
call psb_gefree(x, desc_A,info) call psb_gefree(x, desc_A,info)
call psb_spfree(A, desc_A,info) call psb_spfree(A, desc_A,info)
call mld_precfree(P,info) call P%free(info)
call psb_cdfree(desc_A,info) call psb_cdfree(desc_A,info)
call psb_exit(ictxt) call psb_exit(ictxt)
stop stop

@ -37,22 +37,21 @@
! POSSIBILITY OF SUCH DAMAGE. ! POSSIBILITY OF SUCH DAMAGE.
! !
! !
! File: mld_zexample_ml.f90 ! File: mld_zexample_1lev.f90
! !
! This sample program solves a linear system by using BiCGStab preconditioned by ! This sample program solves a linear system by using BiCGStab preconditioned by
! RAS with overlap 2 and ILU(0) on the local blocks, as explained in Section 6.1 ! RAS with overlap 2 and ILU(0) on the local blocks, as explained in Section 5.1
! of the MLD2P4 User's and Reference Guide. ! of the MLD2P4 User's and Reference Guide.
! !
! The matrix and the rhs are read from files (if an rhs is not available, the ! The matrix and the rhs are read from files (if an rhs is not available, the
! unit rhs is set). ! unit rhs is set).
! !
program mld_zexample_ml program mld_zexample_1lev
use psb_base_mod use psb_base_mod
use mld_prec_mod use mld_prec_mod
use psb_krylov_mod use psb_krylov_mod
use psb_util_mod use psb_util_mod
use data_input use data_input
implicit none implicit none
! input parameters ! input parameters
@ -86,9 +85,11 @@ program mld_zexample_ml
integer :: i,info,j,m_problem integer :: i,info,j,m_problem
integer(psb_long_int_k_) :: amatsize, precsize, descsize integer(psb_long_int_k_) :: amatsize, precsize, descsize
integer :: ierr, ircode integer :: ierr, ircode
real(psb_dpk_) :: t1, t2, tprec, resmx, resmxp real(psb_dpk_) :: resmx, resmxp
real(psb_dpk_) :: t1, t2, tprec
character(len=20) :: name character(len=20) :: name
integer, parameter :: iunit=12 integer, parameter :: iunit=12
type(psb_z_vect_type) :: x_col, r_col
! initialize the parallel environment ! initialize the parallel environment
@ -190,22 +191,28 @@ program mld_zexample_ml
write(*,'(" ")') write(*,'(" ")')
end if end if
! set RAS with overlap 2 and ILU(0) on the local blocks
call mld_precinit(P,'AS',info) ! START SETTING PARAMETER
call mld_precset(P,mld_sub_ovr_,2,info)
! set RAS
call P%init('AS',info)
! set number of overlaps
call P%set('SUB_OVR',2,info)
! build the preconditioner ! build the preconditioner
t1 = psb_wtime() t1 = psb_wtime()
call mld_precbld(A,desc_A,P,info) call P%build(A,desc_A,info)
tprec = psb_wtime()-t1 tprec = psb_wtime()-t1
call psb_amx(ictxt, tprec) call psb_amx(ictxt, tprec)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_precbld') call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_precbld')
goto 9999 goto 9999
end if end if
@ -238,8 +245,7 @@ program mld_zexample_ml
call psb_sum(ictxt,descsize) call psb_sum(ictxt,descsize)
call psb_sum(ictxt,precsize) call psb_sum(ictxt,precsize)
call mld_precdescr(P,info) call P%descr(info)
if (iam == psb_root_) then if (iam == psb_root_) then
write(*,'(" ")') write(*,'(" ")')
write(*,'("Matrix: ",A)')mtrx_file write(*,'("Matrix: ",A)')mtrx_file
@ -257,9 +263,9 @@ program mld_zexample_ml
write(*,'("Total memory occupation for PREC : ",i12)')precsize write(*,'("Total memory occupation for PREC : ",i12)')precsize
end if end if
call psb_gather(x_glob,x,desc_a,info,root=psb_root_) call psb_gather(x_glob,x_col,desc_a,info,root=psb_root_)
if (info == psb_success_) & if (info == psb_success_) &
& call psb_gather(r_glob,r,desc_a,info,root=psb_root_) & call psb_gather(r_glob,r_col,desc_a,info,root=psb_root_)
if (info /= psb_success_) goto 9999 if (info /= psb_success_) goto 9999
if (iam == psb_root_) then if (iam == psb_root_) then
write(0,'(" ")') write(0,'(" ")')
@ -283,7 +289,7 @@ program mld_zexample_ml
call psb_gefree(b, desc_A,info) call psb_gefree(b, desc_A,info)
call psb_gefree(x, desc_A,info) call psb_gefree(x, desc_A,info)
call psb_spfree(A, desc_A,info) call psb_spfree(A, desc_A,info)
call mld_precfree(P,info) call P%free(info)
call psb_cdfree(desc_A,info) call psb_cdfree(desc_A,info)
call psb_exit(ictxt) call psb_exit(ictxt)
stop stop
@ -323,4 +329,4 @@ contains
call psb_bcast(ictxt,tol) call psb_bcast(ictxt,tol)
end subroutine get_parms end subroutine get_parms
end program mld_zexample_ml end program mld_zexample_1lev

@ -42,9 +42,19 @@
! This sample program solves a linear system by using BiCGStab coupled with ! This sample program solves a linear system by using BiCGStab coupled with
! one of the following multi-level preconditioner, as explained in Section 6.1 ! one of the following multi-level preconditioner, as explained in Section 6.1
! of the MLD2P4 User's and Reference Guide: ! of the MLD2P4 User's and Reference Guide:
! - choice = 1, default multi-level Schwarz preconditioner (Sec. 6.1, Fig. 2) !
! - choice = 2, hybrid three-level Schwarz preconditioner (Sec. 6.1, Fig. 3) ! - choice = 1, initialize the default multi-level preconditioner solver, i.e.,
! - choice = 3, additive three-level Schwarz preconditioner (Sec. 6.1, Fig. 4) ! V-cycle with basic smoothed aggregation, 1 hybrid forward/backward
! GS sweep as pre/post-smoother and UMFPACK as coarsest-level
! solver(Sec. 5.1, Fig. 2)
!
! - choice = 2, a V-cycle preconditioner with 1 block-Jacobi sweep
! (with ILU(0) on the blocks) as pre- and post-smoother, and 8 block-Jacobi
! sweeps (with ILU(0) on the blocks) as coarsest-level solver(Sec. 5.1, Fig. 3)
!
! - choice = 3, build a W-cycle preconditioner with 2 Gauss-Seidel sweeps as
! post-smoother (and no pre-smoother), a distributed coarsest
! matrix, and MUMPS as coarsest-level solver (Sec. 5.1, Fig. 4)
! !
! The matrix and the rhs are read from files (if an rhs is not available, the ! The matrix and the rhs are read from files (if an rhs is not available, the
! unit rhs is set). ! unit rhs is set).
@ -90,7 +100,8 @@ program mld_zexample_ml
integer :: i,info,j,m_problem integer :: i,info,j,m_problem
integer(psb_long_int_k_) :: amatsize, precsize, descsize integer(psb_long_int_k_) :: amatsize, precsize, descsize
integer :: ierr, ircode integer :: ierr, ircode
real(psb_dpk_) :: t1, t2, tprec, resmx, resmxp real(psb_dpk_) :: resmx, resmxp
real(psb_dpk_) :: t1, t2, tprec
character(len=20) :: name character(len=20) :: name
integer, parameter :: iunit=12 integer, parameter :: iunit=12
@ -198,36 +209,37 @@ program mld_zexample_ml
case(1) case(1)
! initialize the default multi-level preconditioner, i.e. hybrid ! initialize the default multi-level preconditioner, i.e. V-cycle
! Schwarz, using RAS (with overlap 1 and ILU(0) on the blocks) ! with basic smoothed aggregation, 1 hybrid forward/backward
! as post-smoother and 4 block-Jacobi sweeps (with UMFPACK LU ! GS sweep as pre/post-smoother and UMFPACK as coarsest-level
! on the blocks) as distributed coarse-level solver ! solver
call mld_precinit(P,'ML',info) call P%init('ML',info)
case(2) case(2)
! set a three-level hybrid Schwarz preconditioner, which uses ! initialize a V-cycle preconditioner with 1 block-Jacobi sweep (with
! block Jacobi (with ILU(0) on the blocks) as post-smoother, ! ILU(0) on the blocks) as pre- and post-smoother, and 8 block-Jacobi
! a coarsest matrix replicated on the processors, and the ! sweeps (with ILU(0) on the blocks) as coarsest-level solver
! LU factorization from UMFPACK as coarse-level solver
call mld_precinit(P,'ML',info,nlev=3) call P%init('ML',info)
call mld_precset(P,mld_smoother_type_,'BJAC',info) call P%set('SMOOTHER_TYPE','BJAC',info)
call mld_precset(P,mld_coarse_mat_,'REPL',info) call P%set('COARSE_SOLVE','BJAC',info)
call mld_precset(P,mld_coarse_solve_,'UMF',info) call P%set('COARSE_SWEEPS',8,info)
case(3) case(3)
! set a three-level additive Schwarz preconditioner, which uses ! initialize a W-cycle preconditioner with 2 Gauss-Seidel sweeps as
! RAS (with overlap 1 and ILU(0) on the blocks) as pre- and ! post-smoother (and no pre-smoother), a distributed coarsest
! post-smoother, and 5 block-Jacobi sweeps (with UMFPACK LU ! matrix, and MUMPS as coarsest-level solver
! on the blocks) as distributed coarsest-level solver
call mld_precinit(P,'ML',info,nlev=3) call P%init('ML',info)
call mld_precset(P,mld_ml_type_,'ADD',info) call P%set('ML_TYPE','WCYCLE',info)
call mld_precset(P,mld_smoother_pos_,'TWOSIDE',info) call P%set('SMOOTHER_TYPE','GS',info)
call mld_precset(P,mld_coarse_sweeps_,5,info) call P%set('SMOOTHER_SWEEPS',0,info,pos='PRE')
call P%set('SMOOTHER_SWEEPS',2,info,pos='POST')
call P%set('COARSE_SOLVE','MUMPS',info)
call P%set('COARSE_MAT','DIST',info)
end select end select
@ -236,18 +248,19 @@ program mld_zexample_ml
call psb_barrier(ictxt) call psb_barrier(ictxt)
t1 = psb_wtime() t1 = psb_wtime()
call mld_precbld(A,desc_A,P,info) ! build the preconditioner
call P%hierarchy_build(A,desc_A,info)
call P%smoothers_build(A,desc_A,info)
tprec = psb_wtime()-t1 tprec = psb_wtime()-t1
call psb_amx(ictxt, tprec) call psb_amx(ictxt, tprec)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_precbld') call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_precbld')
goto 9999 goto 9999
end if end if
! set the initial guess ! set the initial guess
call psb_geall(x,desc_A,info) call psb_geall(x,desc_A,info)
call x%zero() call x%zero()
call psb_geasb(x,desc_A,info) call psb_geasb(x,desc_A,info)
@ -275,7 +288,7 @@ program mld_zexample_ml
call psb_sum(ictxt,descsize) call psb_sum(ictxt,descsize)
call psb_sum(ictxt,precsize) call psb_sum(ictxt,precsize)
call mld_precdescr(P,info) call P%descr(info)
if (iam == psb_root_) then if (iam == psb_root_) then
write(*,'(" ")') write(*,'(" ")')
@ -320,7 +333,7 @@ program mld_zexample_ml
call psb_gefree(b, desc_A,info) call psb_gefree(b, desc_A,info)
call psb_gefree(x, desc_A,info) call psb_gefree(x, desc_A,info)
call psb_spfree(A, desc_A,info) call psb_spfree(A, desc_A,info)
call mld_precfree(P,info) call P%free(info)
call psb_cdfree(desc_A,info) call psb_cdfree(desc_A,info)
call psb_exit(ictxt) call psb_exit(ictxt)
stop stop

@ -40,12 +40,9 @@
! File: mld_dexample_1lev.f90 ! File: mld_dexample_1lev.f90
! !
! This sample program solves a linear system obtained by discretizing a ! This sample program solves a linear system obtained by discretizing a
! PDE with Dirichlet BCs. The solver is BiCGStab coupled with one of the ! PDE with Dirichlet BCs. The solver is BiCGStab preconditioned by
! following multi-level preconditioner, as explained in Section 6.1 of ! RAS with overlap 2 and ILU(0) on the local blocks, as explained in Section 5.1
! the MLD2P4 User's and Reference Guide: ! of the MLD2P4 User's and Reference Guide.
! - choice = 1, default multi-level Schwarz preconditioner (Sec. 6.1, Fig. 2)
! - choice = 2, hybrid three-level Schwarz preconditioner (Sec. 6.1, Fig. 3)
! - choice = 3, additive three-level Schwarz preconditioner (Sec. 6.1, Fig. 4)
! !
! !
! The PDE is a general second order equation in 3d ! The PDE is a general second order equation in 3d
@ -68,46 +65,46 @@ contains
! functions parametrizing the differential equation ! functions parametrizing the differential equation
! !
function b1(x,y,z) function b1(x,y,z)
use psb_base_mod, only : psb_dpk_ use psb_base_mod, only : psb_dpk_, done
real(psb_dpk_) :: b1 real(psb_dpk_) :: b1
real(psb_dpk_), intent(in) :: x,y,z real(psb_dpk_), intent(in) :: x,y,z
b1=1.d0/sqrt(3.d0) b1=done/sqrt(3.d0)
end function b1 end function b1
function b2(x,y,z) function b2(x,y,z)
use psb_base_mod, only : psb_dpk_ use psb_base_mod, only : psb_dpk_, done
real(psb_dpk_) :: b2 real(psb_dpk_) :: b2
real(psb_dpk_), intent(in) :: x,y,z real(psb_dpk_), intent(in) :: x,y,z
b2=1.d0/sqrt(3.d0) b2=done/sqrt(3.d0)
end function b2 end function b2
function b3(x,y,z) function b3(x,y,z)
use psb_base_mod, only : psb_dpk_ use psb_base_mod, only : psb_dpk_, done
real(psb_dpk_) :: b3 real(psb_dpk_) :: b3
real(psb_dpk_), intent(in) :: x,y,z real(psb_dpk_), intent(in) :: x,y,z
b3=1.d0/sqrt(3.d0) b3=done/sqrt(3.d0)
end function b3 end function b3
function c(x,y,z) function c(x,y,z)
use psb_base_mod, only : psb_dpk_ use psb_base_mod, only : psb_dpk_, done
real(psb_dpk_) :: c real(psb_dpk_) :: c
real(psb_dpk_), intent(in) :: x,y,z real(psb_dpk_), intent(in) :: x,y,z
c=0.d0 c=0.d0
end function c end function c
function a1(x,y,z) function a1(x,y,z)
use psb_base_mod, only : psb_dpk_ use psb_base_mod, only : psb_dpk_, done
real(psb_dpk_) :: a1 real(psb_dpk_) :: a1
real(psb_dpk_), intent(in) :: x,y,z real(psb_dpk_), intent(in) :: x,y,z
a1=1.d0/80 a1=done/80
end function a1 end function a1
function a2(x,y,z) function a2(x,y,z)
use psb_base_mod, only : psb_dpk_ use psb_base_mod, only : psb_dpk_, done
real(psb_dpk_) :: a2 real(psb_dpk_) :: a2
real(psb_dpk_), intent(in) :: x,y,z real(psb_dpk_), intent(in) :: x,y,z
a2=1.d0/80 a2=done/80
end function a2 end function a2
function a3(x,y,z) function a3(x,y,z)
use psb_base_mod, only : psb_dpk_ use psb_base_mod, only : psb_dpk_, done
real(psb_dpk_) :: a3 real(psb_dpk_) :: a3
real(psb_dpk_), intent(in) :: x,y,z real(psb_dpk_), intent(in) :: x,y,z
a3=1.d0/80 a3=done/80
end function a3 end function a3
function g(x,y,z) function g(x,y,z)
use psb_base_mod, only : psb_dpk_, done, dzero use psb_base_mod, only : psb_dpk_, done, dzero
@ -155,7 +152,8 @@ program mld_dexample_1lev
integer :: i,info,j integer :: i,info,j
integer(psb_long_int_k_) :: amatsize, precsize, descsize integer(psb_long_int_k_) :: amatsize, precsize, descsize
integer :: idim, nlev, ierr, ircode integer :: idim, nlev, ierr, ircode
real(psb_dpk_) :: t1, t2, tprec, resmx, resmxp real(psb_dpk_) :: resmx, resmxp
real(psb_dpk_) :: t1, t2, tprec
character(len=5) :: afmt='CSR' character(len=5) :: afmt='CSR'
character(len=20) :: name character(len=20) :: name
@ -202,30 +200,32 @@ program mld_dexample_1lev
if (iam == psb_root_) write(*,'("Overall matrix creation time : ",es12.5)')t2 if (iam == psb_root_) write(*,'("Overall matrix creation time : ",es12.5)')t2
if (iam == psb_root_) write(*,'(" ")') if (iam == psb_root_) write(*,'(" ")')
! set RAS with overlap 2 and ILU(0) on the local blocks ! set RAS
call mld_precinit(P,'AS',info) call P%init('AS',info)
call mld_precset(P,mld_sub_ovr_,2,info)
! set number of overlaps
call P%set('SUB_OVR',2,info)
! build the preconditioner ! build the preconditioner
call psb_barrier(ictxt)
t1 = psb_wtime() t1 = psb_wtime()
call mld_precbld(A,desc_A,P,info) call P%build(A,desc_A,info)
tprec = psb_wtime()-t1 tprec = psb_wtime()-t1
call psb_amx(ictxt, tprec) call psb_amx(ictxt, tprec)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_precbld') call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_precbld')
goto 9999 goto 9999
end if end if
! set the initial guess ! set the initial guess
call psb_geall(x,desc_A,info) call psb_geall(x,desc_A,info)
call x%set(dzero) call x%zero()
call psb_geasb(x,desc_A,info) call psb_geasb(x,desc_A,info)
! solve Ax=b with preconditioned BiCGSTAB ! solve Ax=b with preconditioned BiCGSTAB
@ -253,7 +253,7 @@ program mld_dexample_1lev
call psb_sum(ictxt,descsize) call psb_sum(ictxt,descsize)
call psb_sum(ictxt,precsize) call psb_sum(ictxt,precsize)
call mld_precdescr(P,info) call P%descr(info)
if (iam == psb_root_) then if (iam == psb_root_) then
write(*,'(" ")') write(*,'(" ")')
@ -275,7 +275,7 @@ program mld_dexample_1lev
call psb_gefree(b, desc_A,info) call psb_gefree(b, desc_A,info)
call psb_gefree(x, desc_A,info) call psb_gefree(x, desc_A,info)
call psb_spfree(A, desc_A,info) call psb_spfree(A, desc_A,info)
call mld_precfree(P,info) call P%free(info)
call psb_cdfree(desc_A,info) call psb_cdfree(desc_A,info)
call psb_exit(ictxt) call psb_exit(ictxt)
stop stop

@ -41,11 +41,21 @@
! !
! This sample program solves a linear system obtained by discretizing a ! This sample program solves a linear system obtained by discretizing a
! PDE with Dirichlet BCs. The solver is BiCGStab coupled with one of the ! PDE with Dirichlet BCs. The solver is BiCGStab coupled with one of the
! following multi-level preconditioner, as explained in Section 6.1 of ! following multi-level preconditioner, as explained in Section 5.1 of
! the MLD2P4 User's and Reference Guide: ! the MLD2P4 User's and Reference Guide:
! - choice = 1, default multi-level Schwarz preconditioner (Sec. 6.1, Fig. 2) !
! - choice = 2, hybrid three-level Schwarz preconditioner (Sec. 6.1, Fig. 3) ! - choice = 1, initialize the default multi-level preconditioner solver, i.e.,
! - choice = 3, additive three-level Schwarz preconditioner (Sec. 6.1, Fig. 4) ! V-cycle with basic smoothed aggregation, 1 hybrid forward/backward
! GS sweep as pre/post-smoother and UMFPACK as coarsest-level
! solver(Sec. 5.1, Fig. 2)
!
! - choice = 2, a V-cycle preconditioner with 1 block-Jacobi sweep
! (with ILU(0) on the blocks) as pre- and post-smoother, and 8 block-Jacobi
! sweeps (with ILU(0) on the blocks) as coarsest-level solver(Sec. 5.1, Fig. 3)
!
! - choice = 3, build a W-cycle preconditioner with 2 Gauss-Seidel sweeps as
! post-smoother (and no pre-smoother), a distributed coarsest
! matrix, and MUMPS as coarsest-level solver (Sec. 5.1, Fig. 4)
! !
! The PDE is a general second order equation in 3d ! The PDE is a general second order equation in 3d
! !
@ -73,46 +83,46 @@ contains
! functions parametrizing the differential equation ! functions parametrizing the differential equation
! !
function b1(x,y,z) function b1(x,y,z)
use psb_base_mod, only : psb_dpk_ use psb_base_mod, only : psb_dpk_, done
real(psb_dpk_) :: b1 real(psb_dpk_) :: b1
real(psb_dpk_), intent(in) :: x,y,z real(psb_dpk_), intent(in) :: x,y,z
b1=1.d0/sqrt(3.d0) b1=done/sqrt(3.d0)
end function b1 end function b1
function b2(x,y,z) function b2(x,y,z)
use psb_base_mod, only : psb_dpk_ use psb_base_mod, only : psb_dpk_, done
real(psb_dpk_) :: b2 real(psb_dpk_) :: b2
real(psb_dpk_), intent(in) :: x,y,z real(psb_dpk_), intent(in) :: x,y,z
b2=1.d0/sqrt(3.d0) b2=done/sqrt(3.d0)
end function b2 end function b2
function b3(x,y,z) function b3(x,y,z)
use psb_base_mod, only : psb_dpk_ use psb_base_mod, only : psb_dpk_, done
real(psb_dpk_) :: b3 real(psb_dpk_) :: b3
real(psb_dpk_), intent(in) :: x,y,z real(psb_dpk_), intent(in) :: x,y,z
b3=1.d0/sqrt(3.d0) b3=done/sqrt(3.d0)
end function b3 end function b3
function c(x,y,z) function c(x,y,z)
use psb_base_mod, only : psb_dpk_ use psb_base_mod, only : psb_dpk_, done
real(psb_dpk_) :: c real(psb_dpk_) :: c
real(psb_dpk_), intent(in) :: x,y,z real(psb_dpk_), intent(in) :: x,y,z
c=0.d0 c=0.d0
end function c end function c
function a1(x,y,z) function a1(x,y,z)
use psb_base_mod, only : psb_dpk_ use psb_base_mod, only : psb_dpk_, done
real(psb_dpk_) :: a1 real(psb_dpk_) :: a1
real(psb_dpk_), intent(in) :: x,y,z real(psb_dpk_), intent(in) :: x,y,z
a1=1.d0/80 a1=done/80
end function a1 end function a1
function a2(x,y,z) function a2(x,y,z)
use psb_base_mod, only : psb_dpk_ use psb_base_mod, only : psb_dpk_, done
real(psb_dpk_) :: a2 real(psb_dpk_) :: a2
real(psb_dpk_), intent(in) :: x,y,z real(psb_dpk_), intent(in) :: x,y,z
a2=1.d0/80 a2=done/80
end function a2 end function a2
function a3(x,y,z) function a3(x,y,z)
use psb_base_mod, only : psb_dpk_ use psb_base_mod, only : psb_dpk_, done
real(psb_dpk_) :: a3 real(psb_dpk_) :: a3
real(psb_dpk_), intent(in) :: x,y,z real(psb_dpk_), intent(in) :: x,y,z
a3=1.d0/80 a3=done/80
end function a3 end function a3
function g(x,y,z) function g(x,y,z)
use psb_base_mod, only : psb_dpk_, done, dzero use psb_base_mod, only : psb_dpk_, done, dzero
@ -163,7 +173,8 @@ program mld_dexample_ml
integer :: i,info,j integer :: i,info,j
integer(psb_long_int_k_) :: amatsize, precsize, descsize integer(psb_long_int_k_) :: amatsize, precsize, descsize
integer :: idim, ierr, ircode integer :: idim, ierr, ircode
real(psb_dpk_) :: t1, t2, tprec, resmx, resmxp real(psb_dpk_) :: resmx, resmxp
real(psb_dpk_) :: t1, t2, tprec
character(len=5) :: afmt='CSR' character(len=5) :: afmt='CSR'
character(len=20) :: name character(len=20) :: name
@ -213,72 +224,61 @@ program mld_dexample_ml
select case(choice) select case(choice)
case(1) case(1)
! initialize the default multi-level preconditioner, i.e. hybrid ! initialize the default multi-level preconditioner, i.e. V-cycle
! Schwarz, using RAS (with overlap 1 and ILU(0) on the blocks) ! with basic smoothed aggregation, 1 hybrid forward/backward
! as post-smoother and 4 block-Jacobi sweeps (with UMFPACK LU ! GS sweep as pre/post-smoother and UMFPACK as coarsest-level
! on the blocks) as distributed coarse-level solver ! solver
call mld_precinit(P,'ML',info) call P%init('ML',info)
case(2) case(2)
! set a three-level hybrid Schwarz preconditioner, which uses ! initialize a V-cycle preconditioner with 1 block-Jacobi sweep (with
! block Jacobi (with ILU(0) on the blocks) as post-smoother, ! ILU(0) on the blocks) as pre- and post-smoother, and 8 block-Jacobi
! a coarsest matrix replicated on the processors, and the ! sweeps (with ILU(0) on the blocks) as coarsest-level solver
! LU factorization from UMFPACK as coarse-level solver
call mld_precinit(P,'ML',info,nlev=3) call P%init('ML',info)
call mld_precset(P,mld_smoother_type_,'BJAC',info) call P%set('SMOOTHER_TYPE','BJAC',info)
call mld_precset(P,mld_coarse_mat_,'REPL',info) call P%set('COARSE_SOLVE','BJAC',info)
call mld_precset(P,mld_coarse_solve_,'UMF',info) call P%set('COARSE_SWEEPS',8,info)
case(3) case(3)
! set a three-level additive Schwarz preconditioner, which uses ! initialize a W-cycle preconditioner with 2 Gauss-Seidel sweeps as
! RAS (with overlap 1 and ILU(0) on the blocks) as pre- and ! post-smoother (and no pre-smoother), a distributed coarsest
! post-smoother, and 5 block-Jacobi sweeps (with UMFPACK LU ! matrix, and MUMPS as coarsest-level solver
! on the blocks) as distributed coarsest-level solver
call mld_precinit(P,'ML',info,nlev=3)
call mld_precset(P,mld_ml_type_,'ADD',info)
call mld_precset(P,mld_smoother_pos_,'TWOSIDE',info)
call mld_precset(P,mld_coarse_sweeps_,5,info)
case(4)
! set a three-level hybrid Schwarz preconditioner, which uses call P%init('ML',info)
! block Jacobi (with ILU(0) on the blocks) as post-smoother, call P%set('ML_TYPE','WCYCLE',info)
! a coarsest matrix replicated on the processors, and the call P%set('SMOOTHER_TYPE','GS',info)
! multifrontal solver in MUMPS as coarse-level solver call P%set('SMOOTHER_SWEEPS',0,info,pos='PRE')
call P%set('SMOOTHER_SWEEPS',2,info,pos='POST')
call mld_precinit(P,'ML',info,nlev=3) call P%set('COARSE_SOLVE','MUMPS',info)
call mld_precset(P,mld_smoother_type_,'BJAC',info) call P%set('COARSE_MAT','DIST',info)
call mld_precset(P,mld_coarse_mat_,'REPL',info)
call mld_precset(P,mld_coarse_solve_,'MUMPS',info)
end select end select
! build the preconditioner
call psb_barrier(ictxt) call psb_barrier(ictxt)
t1 = psb_wtime() t1 = psb_wtime()
call mld_precbld(A,desc_A,P,info) ! build the preconditioner
call P%hierarchy_build(A,desc_A,info)
call P%smoothers_build(A,desc_A,info)
tprec = psb_wtime()-t1 tprec = psb_wtime()-t1
call psb_amx(ictxt, tprec) call psb_amx(ictxt, tprec)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_precbld') call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_precbld')
goto 9999 goto 9999
end if end if
! set the solver parameters and the initial guess ! set the solver parameters and the initial guess
call psb_geall(x,desc_A,info) call psb_geall(x,desc_A,info)
call x%set(dzero) call x%zero()
call psb_geasb(x,desc_A,info) call psb_geasb(x,desc_A,info)
! solve Ax=b with preconditioned BiCGSTAB ! solve Ax=b with preconditioned BiCGSTAB
@ -292,7 +292,7 @@ program mld_dexample_ml
call psb_amx(ictxt,t2) call psb_amx(ictxt,t2)
call psb_geall(r,desc_A,info) call psb_geall(r,desc_A,info)
call r%set(dzero) call r%zero()
call psb_geasb(r,desc_A,info) call psb_geasb(r,desc_A,info)
call psb_geaxpby(done,b,dzero,r,desc_A,info) call psb_geaxpby(done,b,dzero,r,desc_A,info)
call psb_spmm(-done,A,x,done,r,desc_A,info) call psb_spmm(-done,A,x,done,r,desc_A,info)
@ -306,7 +306,7 @@ program mld_dexample_ml
call psb_sum(ictxt,descsize) call psb_sum(ictxt,descsize)
call psb_sum(ictxt,precsize) call psb_sum(ictxt,precsize)
call mld_precdescr(P,info) call P%descr(info)
if (iam == psb_root_) then if (iam == psb_root_) then
write(*,'(" ")') write(*,'(" ")')
@ -328,7 +328,7 @@ program mld_dexample_ml
call psb_gefree(b, desc_A,info) call psb_gefree(b, desc_A,info)
call psb_gefree(x, desc_A,info) call psb_gefree(x, desc_A,info)
call psb_spfree(A, desc_A,info) call psb_spfree(A, desc_A,info)
call mld_precfree(P,info) call P%free(info)
call psb_cdfree(desc_A,info) call psb_cdfree(desc_A,info)
call psb_exit(ictxt) call psb_exit(ictxt)
stop stop

@ -40,12 +40,10 @@
! File: mld_sexample_1lev.f90 ! File: mld_sexample_1lev.f90
! !
! This sample program solves a linear system obtained by discretizing a ! This sample program solves a linear system obtained by discretizing a
! PDE with Dirichlet BCs. The solver is BiCGStab coupled with one of the ! PDE with Dirichlet BCs. The solver is BiCGStab preconditioned by
! following multi-level preconditioner, as explained in Section 6.1 of ! RAS with overlap 2 and ILU(0) on the local blocks, as explained in Section 5.1
! the MLD2P4 User's and Reference Guide: ! of the MLD2P4 User's and Reference Guide.
! - choice = 1, default multi-level Schwarz preconditioner (Sec. 6.1, Fig. 2) !
! - choice = 2, hybrid three-level Schwarz preconditioner (Sec. 6.1, Fig. 3)
! - choice = 3, additive three-level Schwarz preconditioner (Sec. 6.1, Fig. 4)
! !
! The PDE is a general second order equation in 3d ! The PDE is a general second order equation in 3d
! !
@ -61,57 +59,52 @@
! !
! Note that if b1=b2=b3=c=0., the PDE is the Laplace equation. ! Note that if b1=b2=b3=c=0., the PDE is the Laplace equation.
! !
! In this sample program the index space of the discretized module dpde_mod
! computational domain is first numbered sequentially in a standard way,
! then the corresponding vector is distributed according to a BLOCK
! data distribution.
!
module spde_mod
contains contains
! !
! functions parametrizing the differential equation ! functions parametrizing the differential equation
! !
function b1(x,y,z) function b1(x,y,z)
use psb_base_mod, only : psb_spk_ use psb_base_mod, only : psb_spk_, sone
real(psb_spk_) :: b1 real(psb_spk_) :: b1
real(psb_spk_), intent(in) :: x,y,z real(psb_spk_), intent(in) :: x,y,z
b1=1.e0/sqrt(3.e0) b1=sone/sqrt(3.d0)
end function b1 end function b1
function b2(x,y,z) function b2(x,y,z)
use psb_base_mod, only : psb_spk_ use psb_base_mod, only : psb_spk_, sone
real(psb_spk_) :: b2 real(psb_spk_) :: b2
real(psb_spk_), intent(in) :: x,y,z real(psb_spk_), intent(in) :: x,y,z
b2=1.e0/sqrt(3.e0) b2=sone/sqrt(3.d0)
end function b2 end function b2
function b3(x,y,z) function b3(x,y,z)
use psb_base_mod, only : psb_spk_ use psb_base_mod, only : psb_spk_, sone
real(psb_spk_) :: b3 real(psb_spk_) :: b3
real(psb_spk_), intent(in) :: x,y,z real(psb_spk_), intent(in) :: x,y,z
b3=1.e0/sqrt(3.e0) b3=sone/sqrt(3.d0)
end function b3 end function b3
function c(x,y,z) function c(x,y,z)
use psb_base_mod, only : psb_spk_ use psb_base_mod, only : psb_spk_, sone
real(psb_spk_) :: c real(psb_spk_) :: c
real(psb_spk_), intent(in) :: x,y,z real(psb_spk_), intent(in) :: x,y,z
c=0.e0 c=0.d0
end function c end function c
function a1(x,y,z) function a1(x,y,z)
use psb_base_mod, only : psb_spk_ use psb_base_mod, only : psb_spk_, sone
real(psb_spk_) :: a1 real(psb_spk_) :: a1
real(psb_spk_), intent(in) :: x,y,z real(psb_spk_), intent(in) :: x,y,z
a1=1.e0/80 a1=sone/80
end function a1 end function a1
function a2(x,y,z) function a2(x,y,z)
use psb_base_mod, only : psb_spk_ use psb_base_mod, only : psb_spk_, sone
real(psb_spk_) :: a2 real(psb_spk_) :: a2
real(psb_spk_), intent(in) :: x,y,z real(psb_spk_), intent(in) :: x,y,z
a2=1.e0/80 a2=sone/80
end function a2 end function a2
function a3(x,y,z) function a3(x,y,z)
use psb_base_mod, only : psb_spk_ use psb_base_mod, only : psb_spk_, sone
real(psb_spk_) :: a3 real(psb_spk_) :: a3
real(psb_spk_), intent(in) :: x,y,z real(psb_spk_), intent(in) :: x,y,z
a3=1.e0/80 a3=sone/80
end function a3 end function a3
function g(x,y,z) function g(x,y,z)
use psb_base_mod, only : psb_spk_, sone, szero use psb_base_mod, only : psb_spk_, sone, szero
@ -124,7 +117,7 @@ contains
g = exp(y**2-z**2) g = exp(y**2-z**2)
end if end if
end function g end function g
end module spde_mod end module dpde_mod
program mld_sexample_1lev program mld_sexample_1lev
use psb_base_mod use psb_base_mod
@ -132,8 +125,7 @@ program mld_sexample_1lev
use psb_krylov_mod use psb_krylov_mod
use psb_util_mod use psb_util_mod
use data_input use data_input
use spde_mod use dpde_mod
use mld_s_mumps_solver
implicit none implicit none
@ -147,12 +139,11 @@ program mld_sexample_1lev
type(mld_sprec_type) :: P type(mld_sprec_type) :: P
! right-hand side, solution and residual vectors ! right-hand side, solution and residual vectors
type(psb_s_vect_type) :: x, b, r type(psb_s_vect_type) :: x, b, r
! solver parameters ! solver parameters
real(psb_spk_) :: tol, err real(psb_spk_) :: tol, err
integer :: itmax, iter, itrace, istop integer :: itmax, iter, itrace, istop
type(mld_s_mumps_solver_type) :: sv
! parallel environment parameters ! parallel environment parameters
integer :: ictxt, iam, np integer :: ictxt, iam, np
@ -161,8 +152,8 @@ program mld_sexample_1lev
integer :: i,info,j integer :: i,info,j
integer(psb_long_int_k_) :: amatsize, precsize, descsize integer(psb_long_int_k_) :: amatsize, precsize, descsize
integer :: idim, nlev, ierr, ircode integer :: idim, nlev, ierr, ircode
real(psb_dpk_) :: t1, t2, tprec
real(psb_spk_) :: resmx, resmxp real(psb_spk_) :: resmx, resmxp
real(psb_dpk_) :: t1, t2, tprec
character(len=5) :: afmt='CSR' character(len=5) :: afmt='CSR'
character(len=20) :: name character(len=20) :: name
@ -209,37 +200,32 @@ program mld_sexample_1lev
if (iam == psb_root_) write(*,'("Overall matrix creation time : ",es12.5)')t2 if (iam == psb_root_) write(*,'("Overall matrix creation time : ",es12.5)')t2
if (iam == psb_root_) write(*,'(" ")') if (iam == psb_root_) write(*,'(" ")')
! set RAS
call P%init('AS',info)
! set MUMPS as solver ! set number of overlaps
call mld_precinit(P,'AS',info)
call mld_precset(P,mld_sub_ovr_,2,info)
call sv%default
call P%set(sv,info)
call mld_precset(P,mld_mumps_print_err_,10,info)
call P%set('SUB_OVR',2,info)
! build the preconditioner ! build the preconditioner
call psb_barrier(ictxt)
t1 = psb_wtime() t1 = psb_wtime()
call mld_precbld(A,desc_A,P,info) call P%build(A,desc_A,info)
tprec = psb_wtime()-t1 tprec = psb_wtime()-t1
call psb_amx(ictxt, tprec) call psb_amx(ictxt, tprec)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_precbld') call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_precbld')
goto 9999 goto 9999
end if end if
! set the initial guess ! set the initial guess
call psb_geall(x,desc_A,info) call psb_geall(x,desc_A,info)
call x%set(szero) call x%zero()
call psb_geasb(x,desc_A,info) call psb_geasb(x,desc_A,info)
! solve Ax=b with preconditioned BiCGSTAB ! solve Ax=b with preconditioned BiCGSTAB
@ -267,7 +253,7 @@ program mld_sexample_1lev
call psb_sum(ictxt,descsize) call psb_sum(ictxt,descsize)
call psb_sum(ictxt,precsize) call psb_sum(ictxt,precsize)
call mld_precdescr(P,info) call P%descr(info)
if (iam == psb_root_) then if (iam == psb_root_) then
write(*,'(" ")') write(*,'(" ")')
@ -289,7 +275,7 @@ program mld_sexample_1lev
call psb_gefree(b, desc_A,info) call psb_gefree(b, desc_A,info)
call psb_gefree(x, desc_A,info) call psb_gefree(x, desc_A,info)
call psb_spfree(A, desc_A,info) call psb_spfree(A, desc_A,info)
call mld_precfree(P,info) call P%free(info)
call psb_cdfree(desc_A,info) call psb_cdfree(desc_A,info)
call psb_exit(ictxt) call psb_exit(ictxt)
stop stop

@ -41,11 +41,21 @@
! !
! This sample program solves a linear system obtained by discretizing a ! This sample program solves a linear system obtained by discretizing a
! PDE with Dirichlet BCs. The solver is BiCGStab coupled with one of the ! PDE with Dirichlet BCs. The solver is BiCGStab coupled with one of the
! following multi-level preconditioner, as explained in Section 6.1 of ! following multi-level preconditioner, as explained in Section 5.1 of
! the MLD2P4 User's and Reference Guide: ! the MLD2P4 User's and Reference Guide:
! - choice = 1, default multi-level Schwarz preconditioner (Sec. 6.1, Fig. 2) !
! - choice = 2, hybrid three-level Schwarz preconditioner (Sec. 6.1, Fig. 3) ! - choice = 1, initialize the default multi-level preconditioner solver, i.e.,
! - choice = 3, additive three-level Schwarz preconditioner (Sec. 6.1, Fig. 4) ! V-cycle with basic smoothed aggregation, 1 hybrid forward/backward
! GS sweep as pre/post-smoother and UMFPACK as coarsest-level
! solver(Sec. 5.1, Fig. 2)
!
! - choice = 2, a V-cycle preconditioner with 1 block-Jacobi sweep
! (with ILU(0) on the blocks) as pre- and post-smoother, and 8 block-Jacobi
! sweeps (with ILU(0) on the blocks) as coarsest-level solver(Sec. 5.1, Fig. 3)
!
! - choice = 3, build a W-cycle preconditioner with 2 Gauss-Seidel sweeps as
! post-smoother (and no pre-smoother), a distributed coarsest
! matrix, and MUMPS as coarsest-level solver (Sec. 5.1, Fig. 4)
! !
! The PDE is a general second order equation in 3d ! The PDE is a general second order equation in 3d
! !
@ -66,52 +76,53 @@
! then the corresponding vector is distributed according to a BLOCK ! then the corresponding vector is distributed according to a BLOCK
! data distribution. ! data distribution.
! !
module spde_mod module dpde_mod
contains contains
! !
! functions parametrizing the differential equation ! functions parametrizing the differential equation
! !
function b1(x,y,z) function b1(x,y,z)
use psb_base_mod, only : psb_spk_ use psb_base_mod, only : psb_spk_, sone
real(psb_spk_) :: b1 real(psb_spk_) :: b1
real(psb_spk_), intent(in) :: x,y,z real(psb_spk_), intent(in) :: x,y,z
b1=1.e0/sqrt(3.e0) b1=sone/sqrt(3.d0)
end function b1 end function b1
function b2(x,y,z) function b2(x,y,z)
use psb_base_mod, only : psb_spk_ use psb_base_mod, only : psb_spk_, sone
real(psb_spk_) :: b2 real(psb_spk_) :: b2
real(psb_spk_), intent(in) :: x,y,z real(psb_spk_), intent(in) :: x,y,z
b2=1.e0/sqrt(3.e0) b2=sone/sqrt(3.d0)
end function b2 end function b2
function b3(x,y,z) function b3(x,y,z)
use psb_base_mod, only : psb_spk_ use psb_base_mod, only : psb_spk_, sone
real(psb_spk_) :: b3 real(psb_spk_) :: b3
real(psb_spk_), intent(in) :: x,y,z real(psb_spk_), intent(in) :: x,y,z
b3=1.e0/sqrt(3.e0) b3=sone/sqrt(3.d0)
end function b3 end function b3
function c(x,y,z) function c(x,y,z)
use psb_base_mod, only : psb_spk_ use psb_base_mod, only : psb_spk_, sone
real(psb_spk_) :: c real(psb_spk_) :: c
real(psb_spk_), intent(in) :: x,y,z real(psb_spk_), intent(in) :: x,y,z
c=0.e0 c=0.d0
end function c end function c
function a1(x,y,z) function a1(x,y,z)
use psb_base_mod, only : psb_spk_ use psb_base_mod, only : psb_spk_, sone
real(psb_spk_) :: a1 real(psb_spk_) :: a1
real(psb_spk_), intent(in) :: x,y,z real(psb_spk_), intent(in) :: x,y,z
a1=1.e0/80 a1=sone/80
end function a1 end function a1
function a2(x,y,z) function a2(x,y,z)
use psb_base_mod, only : psb_spk_ use psb_base_mod, only : psb_spk_, sone
real(psb_spk_) :: a2 real(psb_spk_) :: a2
real(psb_spk_), intent(in) :: x,y,z real(psb_spk_), intent(in) :: x,y,z
a2=1.e0/80 a2=sone/80
end function a2 end function a2
function a3(x,y,z) function a3(x,y,z)
use psb_base_mod, only : psb_spk_ use psb_base_mod, only : psb_spk_, sone
real(psb_spk_) :: a3 real(psb_spk_) :: a3
real(psb_spk_), intent(in) :: x,y,z real(psb_spk_), intent(in) :: x,y,z
a3=1.e0/80 a3=sone/80
end function a3 end function a3
function g(x,y,z) function g(x,y,z)
use psb_base_mod, only : psb_spk_, sone, szero use psb_base_mod, only : psb_spk_, sone, szero
@ -124,7 +135,7 @@ contains
g = exp(y**2-z**2) g = exp(y**2-z**2)
end if end if
end function g end function g
end module spde_mod end module dpde_mod
program mld_sexample_ml program mld_sexample_ml
use psb_base_mod use psb_base_mod
@ -132,7 +143,7 @@ program mld_sexample_ml
use psb_krylov_mod use psb_krylov_mod
use psb_util_mod use psb_util_mod
use data_input use data_input
use spde_mod use dpde_mod
implicit none implicit none
! input parameters ! input parameters
@ -162,8 +173,8 @@ program mld_sexample_ml
integer :: i,info,j integer :: i,info,j
integer(psb_long_int_k_) :: amatsize, precsize, descsize integer(psb_long_int_k_) :: amatsize, precsize, descsize
integer :: idim, ierr, ircode integer :: idim, ierr, ircode
real(psb_dpk_) :: t1, t2, tprec
real(psb_spk_) :: resmx, resmxp real(psb_spk_) :: resmx, resmxp
real(psb_dpk_) :: t1, t2, tprec
character(len=5) :: afmt='CSR' character(len=5) :: afmt='CSR'
character(len=20) :: name character(len=20) :: name
@ -177,6 +188,11 @@ program mld_sexample_ml
call psb_exit(ictxt) call psb_exit(ictxt)
stop stop
endif endif
name='mld_sexample_ml'
if(psb_get_errstatus() /= 0) goto 9999
info=psb_success_
call psb_set_errverbosity(2)
! !
! Hello world ! Hello world
! !
@ -185,11 +201,6 @@ program mld_sexample_ml
write(*,*) 'This is the ',trim(name),' sample program' write(*,*) 'This is the ',trim(name),' sample program'
end if end if
name='mld_sexample_ml'
if(psb_get_errstatus() /= 0) goto 9999
info=psb_success_
call psb_set_errverbosity(2)
! get parameters ! get parameters
call get_parms(ictxt,choice,idim,itmax,tol) call get_parms(ictxt,choice,idim,itmax,tol)
@ -213,72 +224,61 @@ program mld_sexample_ml
select case(choice) select case(choice)
case(1) case(1)
! initialize the default multi-level preconditioner, i.e. hybrid ! initialize the default multi-level preconditioner, i.e. V-cycle
! Schwarz, using RAS (with overlap 1 and ILU(0) on the blocks) ! with basic smoothed aggregation, 1 hybrid forward/backward
! as post-smoother and 4 block-Jacobi sweeps (with UMFPACK LU ! GS sweep as pre/post-smoother and UMFPACK as coarsest-level
! on the blocks) as distributed coarse-level solver ! solver
call mld_precinit(P,'ML',info) call P%init('ML',info)
case(2) case(2)
! set a three-level hybrid Schwarz preconditioner, which uses ! initialize a V-cycle preconditioner with 1 block-Jacobi sweep (with
! block Jacobi (with ILU(0) on the blocks) as post-smoother, ! ILU(0) on the blocks) as pre- and post-smoother, and 8 block-Jacobi
! a coarsest matrix replicated on the processors, and the ! sweeps (with ILU(0) on the blocks) as coarsest-level solver
! LU factorization from UMFPACK as coarse-level solver
call mld_precinit(P,'ML',info,nlev=3) call P%init('ML',info)
call mld_precset(P,mld_smoother_type_,'BJAC',info) call P%set('SMOOTHER_TYPE','BJAC',info)
call mld_precset(P,mld_coarse_mat_,'REPL',info) call P%set('COARSE_SOLVE','BJAC',info)
call mld_precset(P,mld_coarse_solve_,'UMF',info) call P%set('COARSE_SWEEPS',8,info)
case(3) case(3)
! set a three-level additive Schwarz preconditioner, which uses ! initialize a W-cycle preconditioner with 2 Gauss-Seidel sweeps as
! RAS (with overlap 1 and ILU(0) on the blocks) as pre- and ! post-smoother (and no pre-smoother), a distributed coarsest
! post-smoother, and 5 block-Jacobi sweeps (with UMFPACK LU ! matrix, and MUMPS as coarsest-level solver
! on the blocks) as distributed coarsest-level solver
call mld_precinit(P,'ML',info,nlev=3)
call mld_precset(P,mld_ml_type_,'ADD',info)
call mld_precset(P,mld_smoother_pos_,'TWOSIDE',info)
call mld_precset(P,mld_coarse_sweeps_,5,info)
case(4)
! set a three-level hybrid Schwarz preconditioner, which uses call P%init('ML',info)
! block Jacobi (with ILU(0) on the blocks) as post-smoother, call P%set('ML_TYPE','WCYCLE',info)
! a coarsest matrix replicated on the processors, and the call P%set('SMOOTHER_TYPE','GS',info)
! multifrontal solver from MUMPS as global coarse-level solver call P%set('SMOOTHER_SWEEPS',0,info,pos='PRE')
call P%set('SMOOTHER_SWEEPS',2,info,pos='POST')
call mld_precinit(P,'ML',info,nlev=3) call P%set('COARSE_SOLVE','MUMPS',info)
call mld_precset(P,mld_smoother_type_,'BJAC',info) call P%set('COARSE_MAT','DIST',info)
call mld_precset(P,mld_coarse_mat_,'DIST',info)
call mld_precset(P,mld_coarse_solve_,'MUMPS',info)
end select end select
! build the preconditioner
call psb_barrier(ictxt) call psb_barrier(ictxt)
t1 = psb_wtime() t1 = psb_wtime()
call mld_precbld(A,desc_A,P,info) ! build the preconditioner
call P%hierarchy_build(A,desc_A,info)
call P%smoothers_build(A,desc_A,info)
tprec = psb_wtime()-t1 tprec = psb_wtime()-t1
call psb_amx(ictxt, tprec) call psb_amx(ictxt, tprec)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_precbld') call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_precbld')
goto 9999 goto 9999
end if end if
! set the solver parameters and the initial guess ! set the solver parameters and the initial guess
call psb_geall(x,desc_A,info) call psb_geall(x,desc_A,info)
call x%set(szero) call x%zero()
call psb_geasb(x,desc_A,info) call psb_geasb(x,desc_A,info)
! solve Ax=b with preconditioned BiCGSTAB ! solve Ax=b with preconditioned BiCGSTAB
@ -292,7 +292,7 @@ program mld_sexample_ml
call psb_amx(ictxt,t2) call psb_amx(ictxt,t2)
call psb_geall(r,desc_A,info) call psb_geall(r,desc_A,info)
call r%set(szero) call r%zero()
call psb_geasb(r,desc_A,info) call psb_geasb(r,desc_A,info)
call psb_geaxpby(sone,b,szero,r,desc_A,info) call psb_geaxpby(sone,b,szero,r,desc_A,info)
call psb_spmm(-sone,A,x,sone,r,desc_A,info) call psb_spmm(-sone,A,x,sone,r,desc_A,info)
@ -306,7 +306,7 @@ program mld_sexample_ml
call psb_sum(ictxt,descsize) call psb_sum(ictxt,descsize)
call psb_sum(ictxt,precsize) call psb_sum(ictxt,precsize)
call mld_precdescr(P,info) call P%descr(info)
if (iam == psb_root_) then if (iam == psb_root_) then
write(*,'(" ")') write(*,'(" ")')
@ -328,7 +328,7 @@ program mld_sexample_ml
call psb_gefree(b, desc_A,info) call psb_gefree(b, desc_A,info)
call psb_gefree(x, desc_A,info) call psb_gefree(x, desc_A,info)
call psb_spfree(A, desc_A,info) call psb_spfree(A, desc_A,info)
call mld_precfree(P,info) call P%free(info)
call psb_cdfree(desc_A,info) call psb_cdfree(desc_A,info)
call psb_exit(ictxt) call psb_exit(ictxt)
stop stop
@ -365,4 +365,5 @@ contains
call psb_bcast(ictxt,tol) call psb_bcast(ictxt,tol)
end subroutine get_parms end subroutine get_parms
end program mld_sexample_ml end program mld_sexample_ml

Loading…
Cancel
Save