|
|
@ -563,6 +563,14 @@ program psb_d_pde2d
|
|
|
|
integer(psb_epk_) :: amatsize, precsize, descsize, d2size
|
|
|
|
integer(psb_epk_) :: amatsize, precsize, descsize, d2size
|
|
|
|
real(psb_dpk_) :: err, eps
|
|
|
|
real(psb_dpk_) :: err, eps
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
! Parameters for solvers in Block-Jacobi preconditioner
|
|
|
|
|
|
|
|
type ainvparms
|
|
|
|
|
|
|
|
character(len=12) :: alg, orth_alg, ilu_alg, ilut_scale
|
|
|
|
|
|
|
|
integer(psb_ipk_) :: fill, inv_fill
|
|
|
|
|
|
|
|
real(psb_dpk_) :: thresh, inv_thresh
|
|
|
|
|
|
|
|
end type ainvparms
|
|
|
|
|
|
|
|
type(ainvparms) :: parms
|
|
|
|
|
|
|
|
|
|
|
|
! other variables
|
|
|
|
! other variables
|
|
|
|
integer(psb_ipk_) :: info, i
|
|
|
|
integer(psb_ipk_) :: info, i
|
|
|
|
character(len=20) :: name,ch_err
|
|
|
|
character(len=20) :: name,ch_err
|
|
|
@ -592,7 +600,7 @@ program psb_d_pde2d
|
|
|
|
!
|
|
|
|
!
|
|
|
|
! get parameters
|
|
|
|
! get parameters
|
|
|
|
!
|
|
|
|
!
|
|
|
|
call get_parms(ictxt,kmethd,ptype,afmt,idim,istopc,itmax,itrace,irst,ipart)
|
|
|
|
call get_parms(ictxt,kmethd,ptype,afmt,idim,istopc,itmax,itrace,irst,ipart,parms)
|
|
|
|
|
|
|
|
|
|
|
|
!
|
|
|
|
!
|
|
|
|
! allocate and fill in the coefficient matrix, rhs and initial guess
|
|
|
|
! allocate and fill in the coefficient matrix, rhs and initial guess
|
|
|
@ -615,6 +623,25 @@ program psb_d_pde2d
|
|
|
|
!
|
|
|
|
!
|
|
|
|
if(iam == psb_root_) write(psb_out_unit,'("Setting preconditioner to : ",a)')ptype
|
|
|
|
if(iam == psb_root_) write(psb_out_unit,'("Setting preconditioner to : ",a)')ptype
|
|
|
|
call prec%init(ictxt,ptype,info)
|
|
|
|
call prec%init(ictxt,ptype,info)
|
|
|
|
|
|
|
|
!
|
|
|
|
|
|
|
|
! Set the options for the BJAC preconditioner
|
|
|
|
|
|
|
|
!
|
|
|
|
|
|
|
|
if (psb_toupper(ptype) == "BJAC") then
|
|
|
|
|
|
|
|
call prec%set('sub_solve', parms%alg, info)
|
|
|
|
|
|
|
|
select case (psb_toupper(parms%alg))
|
|
|
|
|
|
|
|
case ("ILU")
|
|
|
|
|
|
|
|
call prec%set('sub_fillin', parms%fill, info)
|
|
|
|
|
|
|
|
call prec%set('ilu_alg', parms%ilu_alg, info)
|
|
|
|
|
|
|
|
case ("ILUT")
|
|
|
|
|
|
|
|
call prec%set('sub_fillin', parms%fill, info)
|
|
|
|
|
|
|
|
call prec%set('sub_iluthrs', parms%thresh, info)
|
|
|
|
|
|
|
|
call prec%set('ilut_scale', parms%ilut_scale, info)
|
|
|
|
|
|
|
|
case default
|
|
|
|
|
|
|
|
! Do nothing, use default setting in the init routine
|
|
|
|
|
|
|
|
end select
|
|
|
|
|
|
|
|
else
|
|
|
|
|
|
|
|
! nothing to set for NONE or DIAG preconditioner
|
|
|
|
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
|
|
call psb_barrier(ictxt)
|
|
|
|
call psb_barrier(ictxt)
|
|
|
|
t1 = psb_wtime()
|
|
|
|
t1 = psb_wtime()
|
|
|
@ -704,13 +731,14 @@ contains
|
|
|
|
!
|
|
|
|
!
|
|
|
|
! get iteration parameters from standard input
|
|
|
|
! get iteration parameters from standard input
|
|
|
|
!
|
|
|
|
!
|
|
|
|
subroutine get_parms(ictxt,kmethd,ptype,afmt,idim,istopc,itmax,itrace,irst,ipart)
|
|
|
|
subroutine get_parms(ictxt,kmethd,ptype,afmt,idim,istopc,itmax,itrace,irst,ipart,parms)
|
|
|
|
integer(psb_ipk_) :: ictxt
|
|
|
|
integer(psb_ipk_) :: ictxt
|
|
|
|
character(len=*) :: kmethd, ptype, afmt
|
|
|
|
character(len=*) :: kmethd, ptype, afmt
|
|
|
|
integer(psb_ipk_) :: idim, istopc,itmax,itrace,irst,ipart
|
|
|
|
integer(psb_ipk_) :: idim, istopc,itmax,itrace,irst,ipart
|
|
|
|
integer(psb_ipk_) :: np, iam
|
|
|
|
integer(psb_ipk_) :: np, iam
|
|
|
|
integer(psb_ipk_) :: ip, inp_unit
|
|
|
|
integer(psb_ipk_) :: ip, inp_unit
|
|
|
|
character(len=1024) :: filename
|
|
|
|
character(len=1024) :: filename
|
|
|
|
|
|
|
|
type(ainvparms) :: parms
|
|
|
|
|
|
|
|
|
|
|
|
call psb_info(ictxt, iam, np)
|
|
|
|
call psb_info(ictxt, iam, np)
|
|
|
|
|
|
|
|
|
|
|
@ -761,6 +789,25 @@ contains
|
|
|
|
else
|
|
|
|
else
|
|
|
|
irst=1
|
|
|
|
irst=1
|
|
|
|
endif
|
|
|
|
endif
|
|
|
|
|
|
|
|
if (ip >= 9) then
|
|
|
|
|
|
|
|
read(inp_unit,*) parms%alg
|
|
|
|
|
|
|
|
read(inp_unit,*) parms%ilu_alg
|
|
|
|
|
|
|
|
read(inp_unit,*) parms%ilut_scale
|
|
|
|
|
|
|
|
read(inp_unit,*) parms%fill
|
|
|
|
|
|
|
|
read(inp_unit,*) parms%inv_fill
|
|
|
|
|
|
|
|
read(inp_unit,*) parms%thresh
|
|
|
|
|
|
|
|
read(inp_unit,*) parms%inv_thresh
|
|
|
|
|
|
|
|
read(inp_unit,*) parms%orth_alg
|
|
|
|
|
|
|
|
else
|
|
|
|
|
|
|
|
parms%alg = 'ILU' ! Block Solver ILU,ILUT,INVK,AINVT,AORTH
|
|
|
|
|
|
|
|
parms%ilu_alg = 'NONE' ! If ILU : MILU or NONE othewise ignored
|
|
|
|
|
|
|
|
parms%ilut_scale = 'NONE' ! If ILUT: NONE, MAXVAL, DIAG, ARWSUM, ACLSUM, ARCSUM
|
|
|
|
|
|
|
|
parms%fill = 0 ! Level of fill for forward factorization
|
|
|
|
|
|
|
|
parms%inv_fill = 1 ! Level of fill for inverse factorization (only INVK)
|
|
|
|
|
|
|
|
parms%thresh = 1E-1_psb_dpk_ ! Threshold for forward factorization
|
|
|
|
|
|
|
|
parms%inv_thresh = 1E-1_psb_dpk_ ! Threshold for inverse factorization
|
|
|
|
|
|
|
|
parms%orth_alg = 'LLK' ! What orthogonalization algorithm?
|
|
|
|
|
|
|
|
endif
|
|
|
|
|
|
|
|
|
|
|
|
write(psb_out_unit,'("Solving matrix : ell1")')
|
|
|
|
write(psb_out_unit,'("Solving matrix : ell1")')
|
|
|
|
write(psb_out_unit,'("Grid dimensions : ",i5," x ",i5)')idim,idim
|
|
|
|
write(psb_out_unit,'("Grid dimensions : ",i5," x ",i5)')idim,idim
|
|
|
@ -775,6 +822,27 @@ contains
|
|
|
|
write(psb_out_unit,'("Unknown data distrbution, defaulting to 2D")')
|
|
|
|
write(psb_out_unit,'("Unknown data distrbution, defaulting to 2D")')
|
|
|
|
end select
|
|
|
|
end select
|
|
|
|
write(psb_out_unit,'("Preconditioner : ",a)') ptype
|
|
|
|
write(psb_out_unit,'("Preconditioner : ",a)') ptype
|
|
|
|
|
|
|
|
if( psb_toupper(ptype) == "BJAC" ) then
|
|
|
|
|
|
|
|
write(psb_out_unit,'("Block subsolver : ",a)') parms%alg
|
|
|
|
|
|
|
|
select case (psb_toupper(parms%alg))
|
|
|
|
|
|
|
|
case ('ILU')
|
|
|
|
|
|
|
|
write(psb_out_unit,'("Fill in : ",i0)') parms%fill
|
|
|
|
|
|
|
|
write(psb_out_unit,'("MILU : ",a)') parms%ilu_alg
|
|
|
|
|
|
|
|
case ('ILUT')
|
|
|
|
|
|
|
|
write(psb_out_unit,'("Fill in : ",i0)') parms%fill
|
|
|
|
|
|
|
|
write(psb_out_unit,'("Threshold : ",es12.5)') parms%thresh
|
|
|
|
|
|
|
|
write(psb_out_unit,'("Scaling : ",a)') parms%ilut_scale
|
|
|
|
|
|
|
|
case ('INVK')
|
|
|
|
|
|
|
|
write(psb_out_unit,'("Fill in : ",i0)') parms%fill
|
|
|
|
|
|
|
|
write(psb_out_unit,'("Threshold : ",es12.5)') parms%thresh
|
|
|
|
|
|
|
|
write(psb_out_unit,'("Invese Fill in : ",i0)') parms%inv_fill
|
|
|
|
|
|
|
|
write(psb_out_unit,'("Inverse Threshold : ",es12.5)') parms%inv_thresh
|
|
|
|
|
|
|
|
case ('AINVT','AORTH')
|
|
|
|
|
|
|
|
write(psb_out_unit,'("Inverse Threshold : ",es12.5)') parms%inv_thresh
|
|
|
|
|
|
|
|
case default
|
|
|
|
|
|
|
|
write(psb_out_unit,'("Unknown diagonal solver")')
|
|
|
|
|
|
|
|
end select
|
|
|
|
|
|
|
|
end if
|
|
|
|
write(psb_out_unit,'("Iterative method : ",a)') kmethd
|
|
|
|
write(psb_out_unit,'("Iterative method : ",a)') kmethd
|
|
|
|
write(psb_out_unit,'(" ")')
|
|
|
|
write(psb_out_unit,'(" ")')
|
|
|
|
else
|
|
|
|
else
|
|
|
@ -825,5 +893,3 @@ contains
|
|
|
|
end subroutine pr_usage
|
|
|
|
end subroutine pr_usage
|
|
|
|
|
|
|
|
|
|
|
|
end program psb_d_pde2d
|
|
|
|
end program psb_d_pde2d
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|